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Summary. Terzaghi’s conception of the nature of one-dimensional soil consolidation 
[1] is shown to lead to a non-linear differential equation. A dimensional analysis of this 
equation and the boundary conditions of the standard consolidation test [2] gives a 
more general explanation of a well known linear relationship between the total con- 
solidation U(¢) after a time ¢ and ¢'”’. By linearizing the equation in a general manner, 
an expression is obtained for U(t) which includes secondary consolidation terms. Two 
solutions of the linearized equation are obtained; the first for the standard consolidation 
test and the second for consolidation under a boundary load increasing uniformly with 


time. 
Introduction. A medium of clay saturated with water and subjected to a constant 


pressure, continues to contract in volume over a long period of time. Terzaghi, in his 


Erdbaumechanik [3], put forward a diffusion theory to explain the time dependent 
nature of this volume contraction. The following three ideas are the basis of his treat- 
ment and will be taken as the starting point of this paper. 

(1) The medium is imagined as a porous water saturated skeleton of particles of 
negligible compressibility, so that a change in the volume of any region of the skeleton 
equals the volume of water displaced through its boundary. 

(2) Stresses applied to the boundary of the clay produce a hydrostatic head h in 
the pore water and a stress at each point of the skeleton of such a direction and magni- 
tude as to maintain static equilibrium at all instants of time. 

(3) It is finally assumed that spatial variations in h cause the pore water to diffuse 
through the skeleton in accordance with Darcy’s law, which equates the quantity of 
water flowing through any small plane surface o in unit time to the product of the 
area of c, the rate of change of A with distance along the normal to ¢ and a constant K 
called the permeability. 

By considering the problem of the consolidation of clay in one direction only, the 
need for a tensor description of the stress distribution is avoided. The one-dimension 
problem is defined by the assumption that the fluid flow and the motion of the skeleton 
particles are in one direction only and that all relevant physical properties are constant 
at any instant of time on any plane with a normal along this direction. The stress func- 
tion is then the pressure P(z, t) which the skeleton in the neighborhood of the plane z 
can support at time ¢ and, like K, will depend on the physical state of the solid matter 
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of the clay. If the skeleton is a homogeneous aggregate of small non-spherical particles, 
the hydrological state of clay is described by two functions. The first, called the void 
ratio e, and defined as the volume of water saturating unit volume of the porous structure 
in the neighborhood of the point, describes the degree of compaction of the aggregate, 
while the second, a distribution function ¢, describes the orientation of the particles 
with respect to the direction of consolidation. It is evident that @ will depend on the 
previous history of the medium but it will be assumed to begin with that this is such 
that K can be treated as a function of « only. 

By assuming a linear relationship between « and P and a constant value for K, 
Terzaghi derived an equation of the form 


K @P/dx® = dP/dt 
for P. 


Let us, however, derive the general mathematical formulation of the ideas outlined 
above and make approximations only when compelled to by the difficulty of solving 


specific problems. 

The equation of one-dimensional consolidation theory. Let us define z to be the 
volume of the skeleton of the clay bounded by a cylinder S of unit cross-sectional area 
with its axis in the direction of consolidation, and the two planes 0, z fixed relative to the 
skeleton. If L(z, t) is the distance between these two planes at time ¢, 


z 


0L/dz =1+e since L= / (1 + 6 dz. 
0 
The quantity of water Q(t) in this region is given by the expression 


Q= [ e dz (1) 


and hence the rate at which water flows through the planes 0, z into the region is 0Q/dt. 


Now by Darcy’s law, 
20 _[ ah /at) . 
ot 0z/ Oz Jo 


[ de J | K ah J 

=e at eh 
Jo ot 1 + € d2_Jo 
The derivative of this equation with respect to z is 


de _ 2 (_K_ 2), 
ot oz\l+e dz 


and so 


Suppose P,(é) is the pressure on two planes z, , 2, , say, bounding the clay. Since 
consolidation is slow, the total internal pressure is also P,(t). Hence 


P(z, t) + h@, ) = Pod), 
and the derivative of this equation with respect to z, 


0P/dz = —dh/dz 
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inserted in (2) gives the one-dimensional consolidation equation in its final form; 


de 89 ( K_ oF) (3 
ati z@ \L He dz] 3) 

Two more equations connecting K, « and P and a complete set of boundary con- 
ditions are needed to specify a soluble mathematical problem. The complexity of the 
geometry makes the problem of predicting the nature of these relationships a difficult 
one, and little can be said beyond the intuitively obvious statement that K tends to 
zero and P increases when e¢ decreases. 

It is shown in the next section that, provided these two relationships do not ex- 
plicitly involve z or ¢, the system can be solved for the case of the consolidation of a 
semi-infinite medium under a constant boundary pressure. In this case a dimensional 
analysis shows that the volume of water displaced from the skeleton is proportional to 
the square root of the duration of the consolidation. 

The square root law. If the functions K and P are dependent on the void ratio 
e alone, Eq. (3) is invariant under the family of similarity transformations 


zZ— az, tat, a> 0. 


The boundary conditions associated with the problem of the consolidation of a semi- 
infinite medium under a constant boundary pressure P, are: 


e =@ att = Oforallz>O0O 


and 


P = P, atz 0 for allt > 0. 


These conditions are also invariant under each member of the set of transformations 
above. If the functions A(e) and P(e) are such that the boundary conditions uniquely 
determine the solution of (3), the functions ¢(z, 4) and ¢(az, a°t) must be identically 
equal for all a > O in the region of positive values of z and ¢t. This implied that ¢ is a 
function of z2/t'’* only, for by writing 2* for az and ¢* for a’t and observing that «(z*, ¢*) 
is independent of @ for positive z and ¢, since in this region it is identically equal to 
e(z, t), we obtain by differentiation, the equation 


de/Oa = 2 de/d2* + Jat de/Ol* = O. 
This equation, when multiplied by a, gives the linear partial differential equation 
2* de/dz* + 2t* de/Ot* = O 
with the general solution 
e = f(e*/t"”) = f(ze/t'”) 


involving the arbitrary function f. 
Let us write y = 2/t'”’, « = f(y). Since 


de/At = (Af/dy)(dy/dt) = —y/2t°”?-df/dy 


and 
6P/dz = 1/t'?-dP/dy 
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Eq. (3) reduces to the ordinary differential equation 


ydf _ d ( K_ a) 4) 


2 dy dy \l + « dy 
while the boundary conditions become: 
P= P,aty =0 


and e¢ tends to ¢) as y tends to infinity. 

The progress of the consolidation is measured by the volume of water which has 
seeped through each unit of area of the surface z = 0 after the application of the boundary 
pressure P, . The total consolidation U(t) in this problem is the limiting value of 
L(z, 0) — L(z, #) as z tends to infinity and hence 


dU/dt = Lim. —dL/dt 


= Lim. —dQ/dl 


z- 


ola et a ate (=) 2 
ite "TH+ eo \dy/ peo 8 
by a direct substitution from Eq. (2). 
K,/(1 + €) depends only on the value ¢, of e which satisfies the boundary condition 


P(e) = Py , while (dP/dy),-. is a constant found by solving Eq. (4) and evaluating 
dP/de-df/dy at y = 0. The differential equation for U(t) can be easily integrated with 
the boundary condition U = 0 at ¢t = 0 to give the result 
: can ‘ Kk dP " 
U’ = 2At’*, where A is the constant F —_—— ‘ (5) 
1+ € dyJ,-o 


This result is of interest since the semi-infinite problem is closely related to the well 
known consolidation problem with the modified boundary conditions 


e=eq oP = PF, = P(e) att = Offord < 2 < ze, 


and 
P = P, for all t > 0 on the planes z = O and z,. 


During the initial stages of consolidation, in fact till P(z,/2t'’*) differs appreciably from 
P; , this is essentially the superposition of two semi-infinite cases; and hence U(t) is 
given initially by 4A?’”. 

Terzaghi’s linear theory of course, also predicts this result as it corresponds to a 


3). Agreement between his formula for U(t) and the experi- 


very special case of Eq. (3). 
mental consolidation curves over the initial stages of compression has no bearing on the 
validity of the assumptions he has used in linearizing the theory. In many cases, agree- 
ment between his theory and exveriment is in this region of linearity only, and the 
conclusion to be drawn is that in this region K and P are effectively functions of « only. 
The experimental determination of K(e) and P(e). If P is a function of « only, 
there is no great difficulty in determining experimentally the functional relationship 
between them. This is usually done by carrying out a sequence of standard consolidation 
tests [2] and calculating the value of the void ratio at the conclusion of each test. The 
problem of determining the relationship between K and e is a more difficult one. 
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The constant A in Eq. (5) seems at first to be proportional to K(¢,). However, this 
is not so, because (dP/dy),-o is indirectly related to the function K(e) by the differ- 
ential equation (4). This equation is insoluble except for simple forms of the functions 
f and k and these forms do not involve enough parameters to describe reasonably the 
functions. An empirical approach which partly avoids solving Eq. (4) is outlined below. 

Suppose a sequence of standard consolidation tests is carried out on the same sample 
0 < z < z, and the load applied to the boundary is increased by a factor 8 at each step. 
If each experiment is carried through to the completion of the primary consolidation, 
and secondary consolidation is negligible, the boundary conditions for the rth test are; 


P = p’"' P, att = Ofor0 < z < 2, 
P = sp" P, atz = Oand z = zg, fort > 0. 

The boundary conditions for the corresponding semi-infinite problem giving the 

gradient of the (U, ¢'””) line are 
P — 8"' P, as y tends to infinity and P = 6’P, at y = 0. 

A measurement of € at the conclusion of each test gives a set of points (e, , B’Po) 
which map the relationship between P and e. Let us choose the values of the three 
constants C’, C’, m in the empirical relation, 

e= C’P* +C” (6) 
which gives the closest fit to the experimental curve. 

If it is assumed that A/(1 + ©) can be approximated by an equation with two arbi- 
trary constants C, n of the form 

K/(1 + © = CP’, (7) 
the auxiliary equation for determining A, , and so the gradient of the (U, t'”) curve 


for the rth test is 


“eee d (c oF) 
FE CP”) = =—(crS). 
2 dy | ' dy ' dy (8) 
Let P 8 ‘p* and y = B?™""*” y’* so that the equation determining A, 
becomes 
c’ y* — a d ( on a) 
C 2 dy* ~~ dy* ‘ dy* (9) 


and the boundary conditions in terms of y and z’ are 
pP* — P, as y* tends to infinity and P* = BP, at y* = 0. 


This system is independent of r and so dP*/dy* has the same value at y = 0 for each 
value of r. We see from (5) the equation for A, is 


_ | K dP dpP* a 
fAi,p = —_—_ — = — wr = = 
1+edP* dy* dy j,-0 ~ 


dy*/ ,.<0 
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or 


r(n + m + 1) 


log A, = —— ~“s log B + log B, 

where B is a constant independent of r. A graph of log A, against r log 8 should therefore 
be a straight line of gradient (n + m + 1)/2 for the values of r over which the forms 
chosen for « and K/(1 + e) are valid approximations. C must now be found by solving 
Eq. (9) and it seems that this must be done numerically. 

Linear theories of consolidation. The solution of a one-dimensional consolidation 
problem is determined theoretically by the experimental relationships connecting 
e, P and K which give an explicit form to (3), and by the appropriate set of boundary 
conditions. In practice, mathematical difficulties make it necessary to linearize (3) in 
order to obtain analytical solutions for boundary conditions more complex than those 
considered above. This had always been done by assuming K/(1 + e) is constant and 
choosing the linear relationship between P and « which best fits the experimental curve. 

Since there are two functions at our disposal, this linearizing process can be carried 
out in a less restrictive manner. It can readily be seen that to linearize (3), it is necessary 
that K, P and « be related by an equation of the form 


— ? 
oe oP —s de (11) 
1 +e 02 Oz 


where k is independent of «. Equation (3) and this relationship, together give a linear 


equation for e of the form 
De 0 Oe 
«= 2 (, 2). (2 
Ot Oz ( Oz } 


If k is constant or a function of ¢t only, (12) can be put in the dimensionless form, 


de Ove wa 
— se (13) 
OT Or 

where rt = [3 k/a’ dt and x = z/a’, a being a convenient measure of length. 7 will be 


used as a measure of time in all that follows and ‘“‘a’’ will usually be half of the volume 
of the skeleton between the boundary planes and the cylinder S of unit cross-sectional 
area normal to these planes. 

Equation (11) gives an explicit expression for K 


K = -k ro (. + ‘) (14) 
which is fortunately not physically unreasonable. e decreases monotonically to a finite 
limit as P increases and so K, as determined by (14), is positive for positive k and tends 
to zero as e decreases. There is still freedom to choose e as any function of P and ¢, as 
one more relationship between e, K and P is required to specify the consolidation problem 
mathematically. It is interesting to note how previous consolidation theories comply 
with (14). Neglecting e’ and taking « as a linear function of P, and k as a constant, gives 
a constant value for K as in Terzaghi’s theory. The linear theory proposed by D. W. 
Taylor and W. Marchant in 1940 [4] has a constant K and an e — P relationship of the 
form 
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e= « — c(P — P.) +  — c’)(P — Po) exp (—nl). 
Equation (14) is satisfied by choosing k a function of ¢ as defined by the equation 
k(t)[e — (ec — c’) exp (—ypd)] = K 
and a simple integration gives 


€ 


=! log k = 6 + (c’ — a. 


ca 

Standard consolidation tests on many types of clays show that « is a function of ¢ 

as well as P, for, when subjected to the constant boundary pressure P, , « continues to 

decrease after the hydrostatic head h should have vanished. Equation (3) is still valid 

for a time-dependent « — P relationship and so this behavior, usually described as 

secondary consolidation, can be incorporated in a linear theory by using a suitable 

e — P — t relationship and a corresponding permeability K, defined by (14). The func- 
tion e(P, ¢) for many soils, can be adequately described by a form 


e(P, t) = B — B’ log P — B” log (1 + 7) (15) 


involving three constants B, B’ and B”. 
k can be chosen as either that constant which makes the curve 


K = Bk [1 + B — B’ log P — B” log (1 + 7))/P 


fit a known relationship between K and P closely, or that which gives the time scale 
rt = kt/a’, which agrees with an experimental curve. 

The solution of the linear consolidation equation with general boundary conditions. 
Let us consider the one-dimensional consolidation of a homogeneous clay medium 
bounded by plane surfaces at z = a and — a and the cylinder of unit cross-sectional 
area, S, normal to these, under a boundary pressure P,(7), which varies during the ex- 
periment. We will also suppose there is an initial variation in the void ratio along the 
axis of the cylinder. In mathematical terms, we wish to find the solution of Eq. (13) 
subject to the boundary conditions 


€).(x) at r = O for z in the interval — 1 < x < 1 


Il 


€ 


and 


e[P.(r), rT] = g(r) say, ont = +1 forr > 0. 


€ 


The Laplace transform of the solution is readily obtained by conventional methods and 


was found to be 


, cosh gz sinh gx 


ex = (gz + J) owe 2 sinh @ — 1s, 


where 
= [ MeL 2 sinh q(1 — 2’) dz’, 


I, = [ cole = ein ql — 2x’) dz’, 
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 eo(x’) . 
I, = [ <= sinh q(x — x’) dx’, 
q 


70 


[ e “"e(x, 7) dr, I, = [ e“'g(r)dr and q =s. 


Jd 70 


€r 


The total consolidation is by definition L(a, 0) — L(a, r) — L(— a, 0) + L(— a, 7), 
where 


| 


L(z, r) = [ (1 + «€) dz 
and hence 
U(r) =a | (e, — €) dx 


while 


r) 1 
U,= [ e"Udr=a [ (s — «.) dx 
Jo oe. 
(16) 


7 2a) f €o(x’) + €o(— 2’) cosh qu’ pee m tanh a]. 


2s cosh q q 


If «,(x) is a constant ¢, say, plus an odd function, the first term of (16) simplifies and 


U,= 2a( - a) tan .. (17) 
Ss q 
In general ¢,(x) affects the behavior of U near r = 0, especially if the deviation from a 
constant plus an odd function is greatest near x = + or — 1. In the standard con- 
solidation test, this effect will appear as an initial departure from linearity in the (U, t'”’) 
curve and so as a form of pre-primary consolidation. 
The Laplace transform (17) can be readily inverted to give 


U= 2a | [eo — g(r — t)]0,(0 | int) dt, (18) 


where @,( | ¢) is the theta function defined p. 388 of [5]. For large values of 7 this integral 
is asymptotic to 2a(e. — g) while for small values of 7 it behaves like 


2anr~*”” [ leg — g(r — O)t'? at (19) 
or like 4a [€, — g(0)] (7/)'”’ if dg/dr is small. Since g(r) = «{[Po(r), r], g(O) is the initial 
void ratio on the boundaries. 

The standard consolidation test. The pressure on the boundaries in the standard 
consolidation test [2] is constant, P, , say, and so g(r) = e(P, , r) is a function describing 
the secondary consolidation. If, for simplicity, we neglect any initial variation in the 
void ratio, we have the immediate result U behaves like 2a(e, — g) for large values of r, 
which is for the relationship (15) 2aB” log (1 + 7). 

Let us consider in more detail, the consolidation of soils adequately described by 
Eq. (15). We see from Eq. (17) that the Laplace transform of the total consolidation is 
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given by the expression 


b= of — Monica], on 
s s q 


where ec, = B — B’ log Py. 

The first term of (20) is the transform of the well known solution of the standard 
consolidation problem, while the second gives the effect of secondary consolidation on 
U, . These terms can be inverted by well known techniques of the Operational Calculus 
to give series of asymptotic solutions useful for small values of + and others for large r. 
Let us write U as U, + U, , where U, is the primary and U, the secondary consolidation 
term. 

For small values of + 


U, = 4a(e. — er M202 12 YS (= 1) [9 77"? exp (—r?/7) — r Erfe o/7)} (21) 


and 
U, = 8aB’n'7§(1 + 7)'? log [7'? + (1 + 9)'?] — 7'7} + 40B” > (-—1)’J, , (22) 
where 


at 


J, = | wt’? exp (—r’/t) log (1 + +r — 2 dt. 
As J, is not expressible in finite terms, the series (22) is not very useful for values of +r 
for which J, cannot be neglected. However, we may use the first term till 

2 log (1 + 7) {x '?7'” exp (—1/7) — erfe (1/7'”’)} 
which is an upper bound to J, , ceases to be negligible. For large values of the dimension- 


less parameter 7 


U, g 2al€, = ah! o-~ = Zz. exp {—(2r a 1)*n’ 7/4], (2r + We (23) 
and 
™=e 2a lg (1 + 7) + 20M, exp [—2°(2r + 1)*2/4] + N,/(1 + ah, 
1 = Q -—2 x a 2(Or 2 Ei* 2 2, 
M x exp [—a (2r + 1)°/4]E1* [x (2r + 1)°/4], (24) 


N, = (-4)°77(4"*? — 1) r Bo 4a/(2r + 4)!, 
where B,,., are Bernoulli numbers and £7*(z) is the logarithmic integral defined p. 2 


of [6] as li (e”). 
The first few terms of the asymptotic expression (24), which is useful for values of 


7 outside the range of (22), are 


, u : 
a 2aB log (1+ 7)- 3(1 +7) i“ 15(1 + 7)? a 


+ .4755 exp (—2.4674r) + 18.895 exp (—22.2077) + -- }. 
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Consolidation under a linearly increasing boundary pressure. It was pointed out 
earlier, that the standard consolidation test will not effectively justify approximations 
which are made to linearize Eq. (3). It is evident that in an experiment designed to test 
an approximation arising from (14), the pressure applied to the boundary planes must 
vary with time. Let us consider the simplest problem of this type, that with a linearly 


increasing boundary pressure. 
Suppose a pressure P, is first applied to the boundaries and maintained constant 
until the void ratio attains a value e, throughout the medium. Now, at r = 0 let us 


apply the boundary pressure 
P(r) = P.(1 + br). 


We see from (15) and (17) that 


Uz = —2a[B’ exp (6/1)Ei(—s/0) + B” exp (@)Ei(—9)] S04 (25) 
= (U3), + (U2), say 


U, is given directly by (22) and (24) and, by modifications in the derivation of these 
formulae, it can be shown that, for small values of 7, 


co 


U; = 8aB'(br)/7{(1 + br)” log [(br)'”? + (1 + bx)'”7]} + 4aB’ D> (—1)'J!, (26) 


r=] 
where 


Ji’=_," | t’’* exp (—r’/t) log [1 + b(r — d] dt 
0 


and .// is bounded above by the expression 
2 log (1 + br)[9™'’’r'” exp (—1/7) — erfe (1/7'”’)], 


while for large values of r 


= 2aB’ = 7. Se Li 
U, ~ 2aB" log (1 + br) + >» M} exp [—(2r + 1)'a'r/4] + N,b'/(1 + br)", 


b r=0 
where JN, is as defined in (24) and 
M! = 89’ exp [—2°(2r + 1)°/4b]Ei* [x°(2r + 1)*/40]. (27) 


It will be found that the two approximate forms for U, and U; are unsatisfactory in an 
interval of values of r. It would be useful to know J, and //{ in this interval and they 
could be calculated by a numerical integration of the integrals defining them, but it 
would be better in this case to use the appropriate form of (18). 

The consolidation of media of compressible particles. It was assumed, in deriving 
Eq. (3), that the particles composing the skeleton of the clay had a negligible com- 
pressibility. This seems to offer one direction in which the theory presented in this 
paper can be generalized, but it is not difficult to introduce a compressibility factor 
into the treatment and show that the conclusions and formulae come out practically 
unchanged. If z is defined as the mass of the skeleton of the clay bounded by the cylinder 
S and the two planes 0 and z fixed relative to the skeleton and ¢ as the volume of water 











1960] ONE-DIMENSIONAL SOIL CONSOLIDATION 347 


saturating unit mass of the porous structure, Eq. (1) is unchanged. If L(z, ¢) is now 
defined by the equation 


. és [ (eo +0 a, 


where 7 is the volume of unit mass of the skeleton, Eq. (2) is also unchanged and Eq. 
(3) is replaced by 
9 4 K oP 
oe _ 2 (_K 2). (28) 
ot dz \n + € 02 
The square root law is now valid if K, P and 7 can be expressed as functions of € 
only. Finally, it can be seen that the linearizing approximation, (11), now becomes 
Kk oP de 
~ ‘ - = —k — 
n + € 02 Oz (29) 
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Group theory and its application to the quantum mechanics of atomic spectra. By Kugene 
P. Wigner. Academic Press, New York and London, 1959. xi + 372 pp. $8.80. 


This is an English translation of the original well known German edition (1931) together with the 
addition of three new chapters. The book still deals entirely with the quantum mechanics of atomic 
spectra. The translator points out that the translation was motivated by the lack of a good English 
work on the subject of group theory from “the physicist’s point of view.” 

There is a particularly striking comment by the author in his preface to the effect that von Laue 
remarked that the most important result of the original edition was the “recognition that almost all 
rules of spectroscopy follow from the symmetry of the problem.”’ 

While this book is written for the mathematical and theoretical physicist it has sections (such as 
chapter 17) which are written in a descriptive fashion for the purpose of showing the reader how to use 
results of group theory without becoming any more involved than necessary with the mathematical 
details. Such sections will certainly appeal to anyone not thoroughly familiar with the subject. At the 
same time the book should appeal to the advanced student of the subject particularly in view of the 
recent additions to the work—additions such as the time inversion transformation. The transformation 
t — —t provides an additional symmetry element which transforms a given state into one where all 
velocities (and spins) are reversed in direction. 

The translation seems to have been very worth while. It should make a considerable difference in 
the number of physicists in the country who learn something about group theory. 

{0HN TRUELL 


Sampling inspection tables. Second edition. By Harold F. Dodge and Harry G. Romig. 
John Wiley & Sons, Inc., New York, and Chapman & Hall, Ltd., London, 1959. 
xl + 224 pp. $8.00. 


In this handsome, enlarged, edition of their much used tables the authors have made many new 
and useful additions. Perhaps the most important is the inclusion of “Operating Characteristic (OC) 
Curves’’. These curves give the probability of accepting a lot, when using a chosen sampling scheme, as 
a function of the percentage of defectives in the lot. Such information is useful in choosing a scheme as 
it is often important to know how a scheme behaves for a range of percentage defectives. 

The introduction to the tables has been expanded so that the procedure and theory of setting up 
single or double sampling schemes are clearly explained. Two general types of consumer protection are 
considered: that based on Average Outgoing Quality Limit (AOQL) and that based on Lot Tolerance 
Precent Defective (LTPD). AOQL schemes ensure that the average quality delivered will, almost 
invariably, be above some chosen limit. LTPD schemes limit the percent of defectives in each sampled 
lot. All those concerned with quality control should find this new edition invaluable. 


W. FREIBERGER 


An introduction to advanced dynamics. By 8S. W. McCuskey. Addison-Wesley Publishing 
Co., Inc., Reading, Mass., 1959. viii + 263 pp. $8.50. 


In the author’s words ‘‘The purpose of this introductory book is to familiarize advanced under- 
graduate students in science and mathematics with a few ideas of classical dynamics not ordinarily 
treated in their courses in elementary mechanics.... The emphasis is on the underlying principles 
and a few simple and familiar applications for illustrative purposes only.” 

The end result is quite a readable account of a rather broad range of topics as indicated by the follow- 
ing list of chapter headings: Fundamentals of Newtonian Dynamics, Hamilton’s Principle and Lagrange’s 
Equations, Central Force Motion, Dynamics of a Rigid Body, Oscillatory Motion, Hamilton’s Equations 


(Continued on p. 360) 
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A DERIVATION OF THE BASIC EQUATIONS FOR HYDRODYNAMIC 
LUBRICATION WITH A FLUID HAVING CONSTANT PROPERTIES* 


BY 
H. G. ELROD** 


Consultant, The Franklin Institute 


Abstract. In this paper small parameter techniques are used to derive Reynolds’ 
lubrication equations, and refinements thereof, from the full Navier-Stokes equation. 
An effort has been made to retain rigor in the development comparable to that used in 
present-day boundary-layer developments. 

To derive the differential equations for flow in a curved film of arbitrary thickness 
requires the use of general tensor analysis. The mathematical manipulations are some- 
what involved, but one of the results—a refined Reynolds equation—can be simply 
written for a journal or slipper bearing as follows: 


af h) opi, 3 ( i) & = Gu 2 ( - +) 
ot ie (i 7 +) 2 T 30 {h b+) ef OO a ~ 


Here: 
D = shaft diameter (infinite for a slipper bearing) 
h = film thickness 
p = fluid pressure 
U’ = shaft surface velocity 
x = distance around shaft in direction of rotation 
z = distance parallel to shaft axis 
u = fluid viscosity 


The error of the above differential equation is of the order of (h/L)*, where L is the film 
length in the direction x. 

1. Introduction. The purpose of this paper is to provide a derivation of the basic 
equations for hydrodynamic lubrication with a fluid having constant properties. It 
is expected that analytical techniques similar to those employed here can be adapted 
to the development of equations applicable to fluids having pressure- and temperature- 
dependent properties. 

The derivation given here applies directly to the geometries of both journal and 
slipper bearings. No difficulties are anticipated in the application of similar analysis to 
other geometries, such as those which can occur in thrust bearings, but the writer wished 
to avoid becoming lost in too much generality. 

In the literature, a close approach to the present work appears in an article by 
Wannier [1]. The present work extends that of Wannier in three respects. First, the 
complete Navier-Stokes equations are used as a starting-point (rather than “Stokes” 
equations). Second, the mean film surface can have finite curvature. Third, the present 
approximation procedure is susceptible to improvement without regression. 


*Received June 4, 1958. 
**Also, Professor of Mechanical Engineering, Columbia University. 
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Figure 1 shows the physical situation to be analyzed. A shaft rotates within a bearing, 
which, for the purposes of analysis, is only partial. A circular shaft of radius, 2, rotates 
at angular velocity, w. The small film thickness, h, between this shaft and the bearing 
surface can be an arbitrary function of position. It is completely filled with a fluid 
possessing constant properties, free to flow in and out of the film where it is exposed 
to ambient conditicns at the edge of the bearing. 

















Fic. 1—Schematic diagram of bearing 


The problem to be posed here is that of predicting the velocity and pressure dis- 
tributions within the film. This problem can perhaps best be stated in terms of the 
variables: 


: R y” . R-r 
Smal 7 = =e ce = > ae 7 » 
L*’ E L ? s h(é', €”) ? (1 1) 


vr 


where L is a characteristic dimension of the bearing. The position of a fluid particle is 
known when its coordinates £’, ¢’ and & are specified. Likewise, its motion is given by the 
time derivatives of these quantities, u’. Thus: 


a / 
ve] (1.2) 


Then it is convenient to define the following set of dimensionless ‘‘velocities”’ 


; Lu‘ 

‘ ‘ 

=. (1.3) 
” v 
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Corresponding to Fig. 1 we note that: 


u, = 0; ¢ = —l, 1=1,2,3. 
uu, = LRw : u, = 0, i= 2,3; = 0. 


v 
We can presume that the local film thickness, h, is given by: 
h=hy exp {e(é', ’) } ’ (1.4) 


where the function ¢ and its derivatives are 0(1). The value of h, sets the magnitude, 
so to speak, of ihe film thickness throughout the bearing 
Finally, we introduce a dimensionless ‘‘pressure” x such that 





r=—?.. (1.5) 
(i) 


The value of z is specified around the periphery of the film, and may there generally 
be taken as zero. 
The significant feature of lubrication hydrodynamics is the smallness of the dimension- 


less parameter 
«e =h,/L. (1.6) 


Now it is possible to imagine at least two experimental sequences in which this parameter 
becomes progressively smaller, while at the same time the boundary conditions on the 
u, and w remain invariant. In the first of these sequences, L is progressively increased 
with h, constant, while R is varied proportionately to L and w is varied as L~*. Thus 
u, remains constant on the shaft, and x on the bearing periphery is unaltered if the 
external pressure distribution is unaltered. In the second sequence, hy is progressively 
reduced and fluids of equal kinematic viscosity ‘“v’’ are employed which have, neverthe- 
less, progressively lower densities; i.e., p ~ hj . (In the event that the peripheral pressure 
is constant, as is usual, then this constant pressure can be taken as datum, and the 
fluid density need not be changed in this second sequence.) Again, the boundary con- 
ditions on the u; and 7 are invariant. The existence of such hypothetical experimental 
sequences is not necessary to the use of (h)/Z) as a small parameter in a mathematical 
expansion, but it does tend to assure that there will be a range of experimental conditions 
for which the early terms of the expansion will provide an adequate description of the 
observed variables. 

Although the boundary conditions of the problem can be made independent of e¢, 
as shown above, the differential equations for the u, and z in terms of £", ”, #*, do contain 
e in such a manner that a non-singular perturbation problem can be formulated. Thus, 
we hypothesize that we can represent the dependent variables as follows: 


u, = Uo’, 2,8) + Uie',#,&) + &UE',#,&) ete. 
r=PC,P,E + PG’, 2) + PLE, &,&). (1.7) 


When these series are substituted into the complete Navier-Stokes equations, and the 
coefficients of the various powers of ¢e are equated to zero, a sequence of differential 
equations is generated for the U; and P, . The lowest order equations are the original 
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equations of Reynolds. Successively better solutions can be found by solving the higher- 
order differential equaticns in a manner explicitly indicated in this paper. If the series 
in e [Eq. (1.7)] converge, true solutions of the Navier-Stokes equation should result. 
2. The metric tensor. The simplest manner of obtaining the differential equations 
for the and 7 is through the use of tensor calculus. Hence we first compute 


u. 
* 
the covariant components of the metric tensor. They are given by: 


s7 Oy! oy’ ‘ 
Jas = z: ona (2.1) 
i=10& OF 
The intermediate algebraic details required to obtain the g. are recorded in the Appendix. 


The resulting matrix is given below 
5 





(1 3 a) (= sy (: y Oh ah i On 
~ R ' AL €' L} at! a& — xe 
Li = Go (E) 2 I 1+ (E 9) ‘25 (2.2) 
L L/ o0& 0é L 0é L” 0& 
h ah ws h ah (ty 
ae dg" . or L 
The determinant of the matrix | gag | 1s: 
F , 3 h\° 472 a’ 
g=|gae| =\1—E ys] Lr. (2.3) 
ki 
The contravariant components of the metric tensor are found from the formula 
an Bites ; 
g° = —-Cofactor of gga in g. (2.4) 
g 
As shown in the Appendix, the contravariant array is as follows 
.s Olnh 
I > of 
eT 0 cc 
. 2 ne 
a 1—-#- 
R Fi 
ee oe , Olnh p 
L?g** = G® = 0 -? = (2.5) 
. 0& 
3 Olnh (2 oink)! 
=—s agl 4 ; 2 ¢ ~agl 4 2 
6! __ ys almh L* , NY ats (y 2) 
(, 3 4) 0é h (1 ep a) 0g 
\ * ~ = 
Finally, the Euclidean Christoffel symbols are required. They are given by the 
formula: 
oi ° el 996 i OD at ao 
Tisé',€,&) = 439 (dee “Pao _ aes) (2.6) 
0& 0& 0& 
or: 
“a ' iol OG, OGac OGaz 
Pis(', #,€) = 44 (222 += 3 ~ 20.1), (2.7) 
at at at 
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Explicit expressions for the I‘, will not be given; rather we shall determine their orders- 


of-magnitude in terms of e. 
We first note that all G*° are 0(e°), with the exception of G**, which is 0(e~’). Thus: 


G** = Ofexp (—26; 4 In 6}. (2.8) 
Next we note that all derivatives of the Gs are 0(e’), with the exception of the derivatives 
of G,, , which are 0(e). Thus: 
OG a8 
ag 
Substituting Eqs. (2.8) and (2.9) into Eq. (2.7), we obtain: 
ri, = Ofexp [(—2 6; 65 + 2 — 5 55) Ine]} 
+ Ofexp [(—2 6} 65 + 2 — 6) 62) In e]} 
+ Ofexp [(—2 6; 65 + 2 — 6) 63) In e]} o=1,2,3 (2.10) 
By picking the lowest power of ¢ for each I", , we reach the following conclusions. 





= Ofexp (2 — 6), 43) Ine}. (2.9) 








i#3, a¥l, B¥1: Ti, = 0) 

i~#3, a=1, or B=1, or a=B6B=1: Pas = O(6)| | (2.11) 
i=3, a#1, or B#1: ri, = O(e) 

~=3, a=B=1: “2 = OC). J 


3. Reduction of the momentum and mass continuity equations. The Navier-Stokes 
equations for a fluid having constant properties are [2]: 


J ar i ; ) o i a ‘ 
0 = rort| ve ae ee (2s +rr.- rire) | 


; ag ag” ag" ag 


ur(% » Tiou ‘) a | gi* op | 
og" p og 


In terms of the dimensionless velccities u, and the dimensionless pressure, 7, this equa- 








(3.1) 


tion becomes: 





= gn * +7 8+nM- 
Oe" OF dé“ ae” oe” 

+ (2s +155. - ri.rss)us | (3.2) 
ag” 


: ee ae 
— ug (2s ~ ring) —e G . 
or" og 


With 7 + 3, we retain in Eq. (3.2) those terms of lowest power in ¢; i.e., of O(e"*) 
and O(e~'). The resulting, simplified equation is: 


yi ps OU aru; du, 5 Ole 1 
0 = -—G 11 ag + ¢ T +l aE): + 4 2713 ~ + Te ae + ae v3 | 





ae (3.3) 
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Explicit expressions for the G** are given by (2.5). The required I's, in Eq. (3.3) are 
given with sufficient accuracy as: 





yi 2 9G 1/7 ( ") yil h (é 
is = = $G — 2 D = —G = 3.4 
Tis “ oe ’ R R o(¢), (3.4) 
ols = O(€) because = = (0. (3.5) 

0& 0& 
2 IG OG 22 2 0G; d 

3 = 1G*| 2—3 — 8 | = O(2) because 33 _ 0 (3.6 
T'33 5, |2 aE 3E° | O(e) because ae )) 


+ age Gn _ D _1(L) 
si ten aG a0 CORR N\A/ OR \h 


By virtue of these relations, Eq. (3.3) becomes 


L\*| 07u; Lh du} h du : Or 3 OW 
ai (4) —375 — 2G" ——+ —e 1G" —=+6G”" o” —2 (3.8) 
h/ (0&)° Ro® Roe ‘ " @g' "ae pt 0g man 


With 7 = 3, we again seek the two lowest powers of ¢ in Eq. (3.2). It is readily found 


that, correct in terms of 0(e€-*) and O(«~*), 


On (3.9) 


Now we examine the mass-continuity equation. In tensor form, it can be written as: 


‘oe 
2 we} + (3.10) 
og" 

where, of course, 

; ny cad . 
ay ee ae (3.11 
Since 1? vanishes at £; = 0 and £, = — 1, we can integrate Eq. (3.10) to get: 
a pot wh , 9 77 h , ; 
og! | (1 —£ PY de" + ae i (1 —£ i de = 0. (3.12) 


This last equation is exact. 
4, Reynolds equation with first correction terms. In this section we shall obtain 


a Reynolds equation valid to 0(h,/L). From Eqs. (3.8) we obtain: 


Ouy as 8 h du, ( ) On 
—;*, — 2G" - 4 - == =G = 1.1 
a(e’)” R 0€° R a€° G ho/ o€'’ (4-1) 
uy _ hh dug _ gn(h)’ a (4.2) 
a’)? RAE” N\A) ak?’ ” 

where: 
1 
G? = ————. = - 1.3 
(1 2) L+2P Rt: alan 
— FR 


G” = 2, (4.4) 
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Retention of terms through O0(e) in Eqs. (4.1) and (4.2) gives: 


ous ~ h OU, = ( . hy ar 
a(t’)? 3 R ag — l + 2& I ho oe! (4.5) 


and 


OUsy h du, (" ) Or : 
ye p  (N) Or i 
ae)” sR Ok” ho} OF oo 


Because, to the degree of approximation presently being retained, 7 is independent 
of £°, Eqs. (4.5) and (4.6) can be integrated with respect to £°. For u, we obtain: 


| _ (LRe | ag) » 4g Bh] JLRo _ 
ue (4 *'? ate + Bee ~ 36 





(4.7) 
3) 2 LRw a ¢3)3 5a: | 
+e +4 + (&) ig |? 
where 
ae i+ On 
1 = Nn 8e 
and terms of 0(¢’) and higher have been neglected. 
For u, we obtain: 
2 gf fU+H) hfe , © aay 
uy a af 2 + R ‘5 + 4 + 6 ’ (4.8) 


where 


— (I-) On 
™ ho ae” 


The “velocities” u, and uf from the above Eqs. (4.7) and (4.8) can now be sub- 


stituted into Eq. (3.12) and the definite integrations called for in this last equation can 
be carried out explicitly. When terms through 0(e) are retained, the result is: 


o fh | | ) (1 _ h Or _ (c - n) Le | 

0& (; ( 2R/ dé R v (4.9) 
0 h ny ) On _ 

+ Oe (H-)) (i. l + oR as alias 


In more conventional variables, this equation becomes: 


ra} ( h ae 0 ‘ h ) Op 0 { ( h )} 
ine ao p oO 2,3 +e) eT =. a ee , 
Ox \/ (1 ) Ox + az |" (1 + D/ dz Gl Ox a! 3D/) ’ (4.10) 


where D is the shaft diameter, x is distance around the shaft, and z is distance along its 


~— 


AaxXIs. 

In all usual applications, h/D « 1, so that Eq. (4.10) is essentially Reynolds’ Lubri- 
cation Equation with correction terms. It has been derived in this paper from the full 
Navier-Stokes equation by a systematic procedure employing « = ho/Z as a small 
parameter. Now: 
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hbk _ (R)(L), “ss 
D h L D ho/ \D 
When the characteristic length Z is maintained constant and D is taken progressively 
larger, h/D — 0, and the equation for a plane slipper bearing is obtained. It is seen that 
Reynolds’ equation for a plane slipper bearing has an error of 0(e’), rather than 0(e). 
5. Development of solutions to the Navier-Stokes equations. The general approxi- 
mation procedure, of which we have just examined the first few steps, will now be 
outlined. 
Substitution of the series 1.7 into Eq. (3.8) gives the following sets of equations 
which result from equating to zero the coefficients of the successive powers of ‘‘e’’ 


3° ls h \? OP, ae ae - 

saa = \7 —1 + Ai, &,€&), (5.1) 

0(&) ho/ 0 

rU; h\? aP 

0 2) ‘ — a | PP eq 3 5 2) 

ae)” ho} dF sin Bit 
OP, y 1 > a 
ap = Ar-ilé, €,€) (5.3) 
Ss 

In these equations, the functions H,_, , /,-, and K,_,; are known functions of space 


involving velocity and pressure functions of lower order obtained earlier in the approxi- 


mation procedure. For example, H,_, (&', &’, &*) would involve 
Ui 5 re A Wes Can xe. 


k—-1 3 


and their derivatives. 
Now from Eq. (5.3) we have: 





1 42 43 - e ee . 
Pe, FE) = Pe, EO) + | — dé’. (5.4) 
0 oO 
Or: 
oP, OP, 0 ae a - 
ae ee te Lad ££) 5.5) 
0g Oo 


In this equation, and hereafter, we shall abbreviate 


P,(0) = P,(¢’, ?’, 0). 


Substitution of Eq. (5.5) into Eq. (5.1) gives: 
aU; h \* aP,(0) "Te we 
A SE i A lie - (2° 3 5.6 
ae)” (?-) ag! MiG, £5). aol 


This last equation can be integrated twice with respect to ¢*. For k > 0, the boundary 
conditions are: 


Ux’, &, 0) = Urge’, #, —1) = 0. (5.7) 
After integration we have: 
ne h \? aP,(0) (€°)? : 
Ue © of) = ( " =— + N,_,(é', , &). (5.8 
: . e) h, og 2 k-1\6 9 § »&) 7.8) 
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A like equation is obtained for Uj . Thus: 


a ' h ) aP,(0) (#*)? oe 
2 1 2 3 _— oll ee AS 7 - . ( 9 
l AE 7 ,&) (m ag” 2 + 0, i(é = my 5 ) 
The mass-continuity equation, Eq. (3.12), gives the following recursion equation 
involving the U}.. 


hui’, &, &) dé’ + [ hUiG', &, &).= Qa, &). (5.10) 


0 


Ire 


Substitution of the U, and U; from Eqs. (5.8) and (5.9) gives: 

0 [h \* dP,(0) a f{h\ OP,(0) 1 

ae" (;-) eet ae (}-) ee Ri€ , €). 6.11) 
Here we have a partial differential equation for P,(0) as a function of t' and é’. It is the 
“diffusion equation” with a variable “diffusion coefficient” (h/h,)* and a source term, 
R,., . For k > 0 the boundary condition to be imposed at the edge of the fluid film is 
PO) = OQ. 

In principle, Eq. (5.11) ean be solved. Then P, (é’, &, &) can be found from Eq. 
(5.4) and U,, (é, #, &) and U; (&’, #, &) from Eqs. (5.8) and (5.9). All information for re- 
peating the same process for U}., , U?,, and P,,, is now available. Provided that the series 

1.7) converge, solutions are developed for the complete Navier-Stokes equations subject 
to the boundary conditions 


ul = U' = LRw f= 0 
Vv 
x= P, = fd, 
where the edge of the fluid film is described by the parametric equations 
f=); & = FO. 
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Appendix. Calculation of covariant components of metric tensor. The transformation 
inverse to Eqs. (1.1) is: 


y' =rsing = (R — £h) sing = (R — #h) sin (42) : (A.1) 


y = Le’, (A.2) 


y’ = R—rcosg = R — (R — &h) cos (4 ). (A.3) 
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From Eqs. (A.1-A.3): 


OY  : , oh 
<—=—(R h) cose — & - in ¢ 
9¢" R ( ~ vy) ¥ S dé" Ss + 
dy’ .3 Oh n 
me | — —Ft Ss el 
ag —_— = 
oy" : 
a. = —hsing, 
0& 
oy” ay” oy” 
v= 0; a=1 <<, = 0 
0& dé 0g 
3 e 
oy L , oh 
— = = (R — &h) sing £° —> cos ¢, 
of R ( s ¥ + ~ 9g" ¥ 
ay” oh 
—s = COs ¢ 
OF OE - 
ay* } 
3 — nh COS¢. 
oe se 


From the foregoing derivatives the covariant components gg can 
according to Eq. (2.1). For example: 
L 3 oh 


L, 3 ile oh , . 9 > o372\¢3 . 
Iu = {B (R —&hr cos’ y + 3g! sin’ g — 2p (R — ENE pet SIN y COS y 


Ss 


R 


L 3 y noe 3 oh\” > L 3 Oh 
iain «a Se sin? oo aia” da ok ae Le: ee OE me 
| ‘Ba ey ee ae Oe OSE ENE ee 
or: 
L > : ? 3 oh\* 
a = 4G =— h) p + \é ae : 
Next 
; oh seme . ah\? : ' , Ol 
Jo. = ( 2h sin’ g + L’ + (: ar) cos gy = L' + (: <5 
0& 0& 0& 


The other g,g can be found by similar summations. 


(A.A) 


(A.5) 


(A.6) 


(A.7) 


(A.8) 


(A.9) 


(A.10) 


be calculated 


(A.11) 


COS ¢ 


(A.12) 


). (A.13) 


All are combined in the matrix 


appearing as Eq. (2.2). The third-order determinant of this matrix is readily found to be 


rae 
| Jae | _ (1 ~¢ Av. 


(A.14) 


Calculation of contravariant components of metric tensor. Equation (2.4) is used 


to calculate the contravariant components, g*”. For example, g 


( 43 a) 4 (E ar 
33 48 + See oe ~~ 2 L oé' 
33 po ‘ ) 3 fy ( oe Pe ) ty 2 S 
g I L | I g R L | 


8 is found as follows: 
(E) dh ah 
L} a&' ae 


| ey a ah (E ay 
Ie OE" OF” died L a¢] | 


(A.15) 
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or: 
_— 1 3 Alnh\? , 1 (e anh)" 
oo he + _ h\? (: a + L? E ag ° (A.16) 
rae) 
fi 
Likewise, we find that: 
my ¢* ah\? (é\? ah ah | 
oe + Eay OBS 
g? = (—1)™ L| (1 -¢}) Le | | . 7 wee a 
| .s h oh 3 h oh | 
|> L* ae L? ag 
or: 
s _ _& alnh, 
@--L% (A.18) 


The matrix in Eq. (2.5) gives the remaining components. 
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BOOK REVIEWS 


(Continued from p. 348) 


and Phase Space, The Hamilton Jacobi Equation. However, the main purpose of the book, i.e. to place 
the emphasis on the underlying principles, cannot be regarded as fulfilled. In fact the underlying prin- 
ciples receive a rather cursory treatment and the presentation of some of them lacks clarity and 
precision and may be misleading. It may suffice to cite one instance. It is stated on p. 44: “Suppose 
there are n particles in the system and that F; represents the resultant force acting on the jth particle. 
Then the principle of virtual work states that the system will be in equilibrium when 


A 
> F;- 67; = 0. (2-5) 
j7=1 

Here 6r; is the virtual displacement of the jth particle and is arbitrary except for the constraint con- 

ditions.”’ Of course, if F; is indeed the resultant force, then 57; need not be compatible with constraints. 

By F; the author undoubtedly meant the resultant applied force. But nowhere in the book are the 

applied forces distinguished from the forces of constraint. In fact, constraint forces such as string tensions, 


surface supports are listed as applied forces on p. 48 of the book. 
E. T. Onat 


Mathematics in physics and engineering. By J. Irving and N. Mullineux. Academic Press, 
New York and London, 1959. xvii + 883 pp. $11.50. 


The following table of contents indicates the scope of the work: 


. Introduction to partial differential equations (68 pp.). 

. Ordinary differential equations: Frobenius’ and other methods of solution (58 pp.). 
. Bessel and Legendre functions (80 pp.). 

The Laplace and other transforms (45 pp.). 

. Matrices (59 pp.). 

. Analytical methods in classical and wave mechanics (55 pp. ). 

. Calculus of variations (70 pp.). 

. Complex variable theory and conformal transformations (64 pp.). 
. The calculus of residues (82 pp.). 

10. Transform theory (73 pp.). 

11. Numerical methods (65 pp.). 

12. Integral equations (56 pp.). 


OWONMOARWHe 


There is an appendix (84 pp.) of miscellaneous elementary topics in pure mathematics which supple- 
ments the text. 

The work is suitable as a basis for a first year graduate course on mathematical methods for physicists 
and engineers. It is clearly and concisely written, exceptionally well printed and contains a wealth of 
worked and unworked examples with solutions taken from many areas of applied mathematics; the 
emphasis in both text and the examples is not on mathematical nicetie~ but on applications to, for 
instance, elasticity, supersonic flow, electromagnetism, wave mechanics and heat flow. There is a useful 


list. of general references at the end of each chapter. 
WALTER FREIBERGER 


Fachbegriffe der Programmierungstechnik. Edited by J. Heinhold. R. Oldenbourg Verlag, 
Munchen, 1959. 34 pp. $1.05. 


There has been no field in recent years in which questions of terminology have been as confusing 
as in digital computer science. The present booklet is published by the Gesellschaft fiir Angewandte 


(Continued on p. 874) 
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AN IMPLICIT, NUMERICAL METHOD FOR SOLVING THE 
TWO-DIMENSIONAL HEAT EQUATION* 


BY 
GEORGE A. BAKER, Jr. 
AND 


THOMAS A. OLIPHANT 
Los Alamos Scientific Laboratory, University of California, Los Alamos, New Mexico 


1. Introduction. We develop an implicit scheme for the numerical solution of the 
two-dimensional heat-flow problem. In the linear case we are able to solve exactly the 
full two-dimensional set of implicit equations. This solution is possible because we 
choose a difference scheme for which the equations are factorable into two one-dimen- 
sional sets. This factorability is basic to the method. 

We extend this method to non-linear equations and non-rectangular regions by the 
use of an iterative scheme to solve the implicit equations obtained. This scheme provides 
second-order convergence, and in the cases we have tested only a very few iterations 
per time step were required. 

The method is proved to be unconditionally stable both in the linear and non-linear 
cases. (We consider only a special set of non-linear problems in the stability analysis.) 
We prove stability in the linear case by the usual type of Fourier analysis and super- 
position of solutions. In the non-linear case we show that the norms of the solutions of 
the difference equation with homogeneous boundary conditions tend to zero as time tends 
to infinity. This method is limited in mesh size and time-step length only by the require- 
ments of accuracy and not stability. 

By way of illustration we include a discussion of two numerical examples. 

2. The linear case. The basic partial differential equation which we wish to con- 
sider is 

0, 90 _ -2 00 (2.1) 

an? ay * at’ 
where x and y are space coordinates, ¢ is time, and a is a constant. We wish to approxi- 
mate this differential equation by a finite difference scheme to allow the approximate 
numerical calculation of the function @. If we denote spatial points by either the pair 
of indices (zj) or the pair (kl), and time by n, then we may write, using the Einstein 
summation convention, a general, linear difference scheme as 


Bui "Oi; = "Au , (2.2) 


where the "a,, depend on quantities with time earlier than n. We wish to choose a par- 
ticular Bi} which will both represent (2.1) and allow for easy calculation of "6,, for 
successive values n. Clearly, an explicit differencing scheme (B*' a 5:6; , where 6: is the 
Kronecker delta) has both these properties, however, such a scheme is well known to 
have the disadvantage of being only conditionally stable. Implicit schemes are usually 


*Received August 25, 1958. Work performed under the auspices of the U. S. Atomic Energy Com- 


mission. 
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unconditionally stable. One such implicit scheme is the Douglas-Peaceman alternating- 
direction method [1]. It represents the application of the one-dimensional Bruce, Peace- 
man, Rachford, and Rice [2] method to two dimensions. It has the disadvantage that 
only one direction is treated exactly at each advance in time. We feel it desirable to 
choose a differencing scheme which will permit an exact treatment over the entire two 
dimensional mesh. With this view in mind we note that if we may factor B;} as 

nn = A,B; (2.3) 


then we may triangularly resolve A; and B; separately as was done in one dimension 
by Bruce, Peaceman, Rachford, and Rice and obtain a method for the direct calculation 
of the "@,; from ”@,; with m < n. We shall now restrict ourselves to 9-point difference 


schemes, i.e., 

Bii=0 if |t-—k|>1, or |j-—1|]>1. (2.4) 
It will be shown later that 5-point difference schemes are unfactorable. Thus let A; and 
B; be triple-diagonal matrices, 


Bi = (y: 8:°' + yo 61 + ys 837’). 


++1 ' et et-l 
A, = (x, 6; + 2X. 6. +73 & ), (2.5) 


In order to represent Eq. (2.1), B;; must be a numerical representation of the Laplacian 
| ’ l I 
plus whatever part of the representation of the time derivative involves 6 at the advanced 
time. Thus (suppressing the superscript 7), 
Biié;; represents 80. + V7 6x: . (2.6) 
From the symmetry properties of the Laplacian, we may restrict the values of Bj} 
(allowing a different mesh-spacing in the x and y directions) to 


k+1,1 k-1,l 

By = DEI = By 

|? hl = tee = r sg 
via - : (2.7) 
k+1,1+1 k+1,1—1 k—-1,1+1 k-1,1—1 

Bi; “= kL = Bui = By; By 

Bi = £. 


Expanding (2.3) by use of (2.5) we get from (2.7) 


n = LM3Y2 = We , 


N= Yl2 = Yoitr ; (2.8) 
B= QYi = G3Yi = VH1Ys = G3Ys ; 
E = X22 - 
These nine equations in six unknowns imply the restrictions 
mv = né, (2.9) 
Y= 243, Yi = Ys. (2.10) 


It should be noted that if u is zero, as for a five-point scheme, then by (2.9) either \ or 7 
must also be zero so that no factorable five-point scheme is possible. 
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More generally it can be shown that if we let 
mn 


Bi, = wWirlan 


with w"*? = Owhenk —-m>1,k—m<0,l1—n>1,orl-—n< 0 and b,i = 0 
wheni—-m>1,i—m<0,j —n>1,orj —n <0, then either (2.9) or = 4u 


must hold. 
To obtain a proper representation of the Laplacian, we consider the Taylor series 


expansion of 8, 
Ax + Ar, y + Ay) = (x, y) + bAx + cAy + d(Az)’ + e(Az)(Ay) + f(dy)’ (2.11) 
+ g(Az)* + h(Az)? Ay + i(Az)(Ay)’ + j(Ay)*. 


Applying the difference scheme (2.7) to (2.11) we obtain 


Bii0,; = (4u + 2n + 20+ HO. + (4u + 2m) d(Az)? + (4u + 2d)f(Ay)’. (2.12) 


As 
a a = 2d +f), (2.13) 
we must have, by (2.6), 
(Qu + »)(Azx)’? = 1, (2.14a) 
(2u + d)(Ay)? = 1, (2.14b) 
4u +2n +2. +£ =B8. (2.14¢) 
Solving (2.9) and (2.14) we obtain 
u = 1/[B(Ax)*(Ay)’), (2.15a) 
n = (Ax) *{1 — 2/[8(Ay)"}}, (2.15b) 
\ = (Ay)~?{1 — 2/[6(Az)’]}, (2.15¢) 
£ = B{1 — 2/[B(Ay)*}} {1 — 2/[B(Az)*}}. (2.15d) 


We may now solve (2.8) and (2.15) for x; and y,; . Thus, 
ay ras af)?2171§ get) al a i i-1 
be : |B( Az) Ay ] ) O- +. [B(Ax) 2] 6: + Oy } | , (2.16) 
x {6:°' + [B(Ay)? — 2] 6 + 8r"}. 
Or, 
Bii = [B(Ax)*(Ay)?)'AiBi , (2.17) 
where, letting 
P = / 2 oun 
u = B(Az) 2, (2.18) 
v = B(Ay)* — 2, 
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we define fu 1000 
«twee 


0 
001 ei (2.19) 
0 


























and fy 1000 ] 
ls £ 8 0 
014 1 +0 
Mal Ft? ? (2.20) 
08 60 Tt 
We may factor A; and B; as follows: 
Tu — & 0 0 0 0 1 
| u— 0 0 0 
0 1 = 82 0 0 
AS = Wy exe = 0 0 1 u— & 0 
0 0 0 l Uu— 8; 
| * 
= _ (2.21) 
fl s, 0 O O | 
01s 0 
| 001% | 
| 
x 10 0 O0O 1 & 
oct 8. I ? 
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where 
& = 0 8 = (u — &-,)7' (2.22) 
and 
[v — ro 0 0 0 0 | 
1 v—?N 0 0 0 
0 1 v—? 0 0 
i ot ew 0 0 1 v—Ts 0 
0 0 0 l v—T% 
| J 
1 r, 0 ? (2.23) 
6 it 0 
00 1 r O 
x 0001f ; 
00001 
where ; ; 
mo=0 m=@—M%-)". (2.24) 
If we now define 
"Jou = 2b, bi, "8:; (2.25) 
and 
Wer = Wy Wr (2.26) 
then (2.2) becomes, by (2.17), (2.21), and (2.23), 
Wit “Guu = B(Az)*(Ay)’ "an: « (2.27) 


Now as w';; involves only values of (vw) such that v < k, w < 1, we may proceed from 
low index numbers to high ones and solve (2.27) for the "g,, by straightforward elimi- 
nation. Then as ,b; ,b; involves only values of (7j) such that ¢ > v, 7 > w, we may proceed 
from high index numbers to low ones and solve (2.25) for the "6;; by straightforward 
elimination. The formulae involved may be written conveniently as 





n B(Az)*(Ay)? "Q;; — “Qi-r5-1 — (v — T;-1) "Oi-1.; - (u _ 8;-1) "Oii-1 
n= ; : > 2.28) 
g @ — ra) — 8) — 
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na n ae — . n oO¢ 
6; = ii = oa 65 41,; esas f 0; 541 — 67; Oy41, 541 (2.29) 


with appropriate modifications at the boundaries. 
3. Generalization to the non-linear case. In this section we shall consider generaliz- 
ing the method given in Sec. 2 to the solution of non-linear partial differential equations 


of the type 
ay , ay _ ] 


Ow . 
oe oe aw p) - ¥ ‘ (3.5) 
OX OY 


ne. 
Y, ¥ at 
In the linear case we modified the equation under consideration by adding 6y to each 


side. Here we obtain 
2 —7 2 ) dy ‘ \ 
By + V*v = BY + h(z,y, WH = (3.2) 
tf 


In the linear case we chose 8 so that the right-hand side was independent of @ at the 
advanced time. Here this is not possible. Let us, however, choose 


0>6> -r Minimum d iC, 4,, va + v: nsed | hs (3.3) 
ii OW |i: 

where X is the coefficient of y at the advanced time in the representation of the time 
derivative. In the heat flow problem h will be positive and in most cases of interest (see 
Sec. 6) will vary roughly like y’, where 0 < y < 1. Let‘us now guess a value of y, and 
substitute it in the right side of the difference equation approximation to (3.2). We may 
now calculate by the method of Sec. 2 the values of ¥;; which appear on the left side. 
If our guess was close to the true solution, we may expand everything to first order in 


the error. Let 


‘lated — Verue 7 


bo 


Then, by (: 


Il 


' es ah | 7 
Be;; + Vie E + Ah(2j) +X A, se |e ; (3.5) 
€ : 


where \ A;; is the representation of the time derivative. Let us restrict the possible choice 
of representations of the time derivative such that (all y > 0) 


(Ai)/C¥iis) <1. (3.6) 

This condition will be satisfied if, for instance, 
n n n—1 9 - 
Aj; = Vii a Vii; (3.7) 


and it is satisfied for our choice (4.1) as long as the temperature does not drop by a 
factor of more than 4 at one time step. We shall neglect V“e;; as small in (3.5) because 
we expect our guess to have only small systematic errors, rather than the type of errors 
which would create large errors in the second derivative. Thus, 


( 


o 9} : 
6:5 & 1 + | sca + hy rv Violet = —x,;€%, . (3.8) 
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Because of our selection (3.3) of 8, x must be positive or zero. Therefore, if we take a 
weighted average Of Wpuess ANd Peatoutatea 2S 
Po Westiiiadei + ) 
vy ee ’ 
‘+z 
then y’ will agree with ¥,,,. to within second order in e«*. If y’ is now used as a new 


guess, we may continue our iteration procedure and be assured of second-order conver- 
gence. We shall see in Sec. 6 that this scheme provides very rapid convergence in a sample 


(3.9) 





case. 

4. Stability. Let us first consider the stability of our method as described in Sec. 2 
for the case of linear flow. We pick a representation of 06/d¢ which, to within terms of 
third order, gives 00/dt evaluated at the advanced time, namely, 

(1.50" — 20""* + 0.50""”)/At. (4.1) 
From (2.1) and (2.6) we see that 6 of Sec. 2 is 
B = —3/(20’ At). (4.2) 


Let us obtain the exact solutions to our difference equations. We may expand any 
function on the mesh points in a Fourier series, so it will be sufficient to consider the 


behavior of 
"0;; = 0) exp (ib,x; + tbey;)X(n). (4.3) 
We shall assume that 
X(n) = 2*"*, (4.4) 
Substituting (4.3) into (2.1) as expressed by (2.16) we obtain, making use of trigono- 
metric half-angle formulae, 
lr = [—8(Az)? + 4 sin? (4b, Ax)][—8(Az)” + 4 sin’ ($b, Ay)](6 Ax Ay)? (4.5) 
(4Z-* — Z)/3. 
Thus, 
Zz =24(4-3r)”. (4.6) 
By definition, as 6; and 6, are real, f > 1. Thus, for T < 4/3, Z. are both real and 
I1>Z>342>2, >}. (4.7) 
If ! > 4/3, then Z. are both complex and 


Z.| = (31)"” < }. (4.8) 





Therefore as | Z. | < 1, the difference scheme is unconditionally stable. 
In the limit as Az, Ay, and At tend to zero, 


(Z_)‘/4' — exp [—a’°(bi + b3)¢], (4.9) 
(Z,)'/4' +3°'*' exp [4a°(bi + b2)é]. 


The root Z_ represents the analytic solution of (2.1). The root Z, enables us to represent 
any computational error involved in advancing from time n — 2 ton — 1. We see from 
(4.7) and (4.8) that this error is damped by at least a factor of 2 at each time step. 
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The behavior of the solution of the difference equations in the limit as t — © with 
At + 0 is also of interest. This behavior may be obtained by setting Z = 1 in (4.5) and 
removing the reality requirement on b, and b, . Thus, for the steady state 


(Ax)*[4 sin’ (2b, Ax)]~' + (Ay)*[4 sin’ (3b, Ay)]7' = —2a’ At/3. (4.10) 


If we have enough mesh points so that Arb, « 1, Ayb, « 1, then we may expand the 
sines to first order and by some manipulation obtain 


(bi + b5)/bi  Qbia’® At/(3 + 2b,a’ At). (4.11) 
In the differential equation, the steady state solution is characterized by 
b? + b? = 0. (4.12) 


In order to get a good steady solution we see from (4.11) and (4.12) that bja° At must 
be small. In order to get a good time-dependent solution we must be near the limit given 
in (4.9). This requires that b, Az and b, Ay be small, and bja’At be small. Hence, we will 
get a good asympotic solution when the requirements are met for a good time-dependent 
solution. We see from the above analysis that this method is limited in mesh size and 
time-step length only by the requirements of accuracy and not by stability. 

In the non-linear case, the analysis is more difficult but proceeds in a similar manner. 
We shall not investigate stability in the general case (3.1) but shall restrict A to be of 
the form 

h(x, y, v) = f(z, yv’," (4.13) 


where 0 > y > — 2. We shall for convenience consider only the case of homogeneous 
boundary conditions. We can, of course, simulate non-homogeneous boundary con- 
ditions by letting f be very large near the boundary and so make a region near the 
boundary an effective heat reservoir. In this analysis we follow Bellman [3] and extend 
his analysis to the non-linear case. We remark that it is both necessary and sufficient 
for stability in the usual sense that 

limit >> AC, j, "¥i)C¥ii)” = 0 (4.14) 


hold for any initial conditions and homogeneous boundary conditions. We shall prove 


that our differencing scheme is unconditionally stable in the sense that (4.14) holds for 


any [At/(Az)*]. 

Let us fix our attention on a time, which we shall denote by n. Let us attempt to 
separate variables. From (2.17) and (3.1) we obtain the equation which must be satisfied 
for some separable component, v. It is 
{ (B(Ax)?(Ay)?]"'AiBi — B 4; 81} "W(t h(k, l, Wi (v) 4.15) 


where we assume that 
(4.16) 


<< 
N 
<= 


and that 


"vi = z. "yi; (v)p(v). (4.17) 
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In fixing our attention on a time n we imagine that "y,; is somehow known. If we set 
[3 — 4Z-"(v) + Z-*(v)]/(2 Ad) = dr (4.18) 


then it is easy to show that the solution of (4.15) is the same as the solution of the follow- 
ing problem. Let (suppressing the n and v) 


, yi oF Bimeg Wer , West — We P 
ra F [gs + (age 
S ( Aa Ay (4.19) 
— (Wrsr rer = Weer — Wearer + var)’ 


B(Ax)*(Ay)” 


and 


G = D> hk, l, ved) (ves)? = 1. (4.20) 
k,l 


Find the extrema of F + AG subject to (4.20) and homogeneous boundary conditions. 
If we have N X M interior mesh points, then the theory of quadratic forms [4] assures 
us, as 8 is negative, that there exist NM orthogonal vectors ["y,;(v)] which satisfy (4.19) 
and (4.20) and therefore (4.15). The orthogonality is in the sense that 


("Werlv), "Wer(u)) = 2 h(k, 1, "Wer) "Werlv) "Yar(u) = 6, - (4.21) 


As \ can be shown to be the negative of 7, we must have 
<0. (4.22) 


We may now solve (4.18) for the corresponding values of Z.(v). Equation (4.18) implies 
that for the T of (4.5) we get 


r=1-—-2At/3>1 (4.23) 


by (4.22). Thus | Z. | < 1 as in the linear case. Let us now expand “y,,; in terms of our 
set of orthogonal vectors ["y,;(v)]. In order to reproduce ""'y;; and "~“y,, it will be neces- 
sary to use parts proportional to Z, and to Z_ ; however, we have enough freedom to 
effect this decomposition as "y;; is assumed to be calculated by our difference scheme 
from the values at the two previous times. Let us now compute the time derivative of 


the Norm of "y,; 


Sf kt, Wd Wid" = DICK, D = (Was) 
Cc 


ot jl 


(4.24) 





0" 
= (2 a Y) y h(k, l, "Wer “War i te 
k,l 


If we now approximate dy/dt by an expression of the form of (4.1) and use our ex- 
pansion in terms of "y,;(v), then, if p(v) denotes the coefficient of "y,;(v) in "Y,; , we get 


* Norm ("¥;;) = (2+ 7) p> | p(v) |* AC) (4.25) 
< Max [A(v)]}(2 + 7) Norm ("y,;;). 


c 


We must now examine the behavior of Max, [A(v)]. First let us ecnsider the maximum 
attained by Max, [(A(v)] with respect to all variations of components for a fixed Norm. 
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As this represents continuous variation over a closed and bounded set, the maximum of 
Max, [A(v)] must be attained at some point of the set by Weierstrass’ theorem. However, 
applying the arguments that lead to (4.22) we see that for fixed Norm, 


A(v) < A(Norm) < 0. (4.26) 


We may consider a variation in Norm by simply scaling all the y,,; . It easily follows 
from (4.19) and (4.20) that for any Norm 


A(v) < Ao [Norm ("y;,;)]~”. (4.27) 
Thence, by (4.25) and (4.27) we have 
4] . . J, 3 " Ss 
~ [Norm ("y;;)] < Ao(2 + y) [Norm ("y;,;)]'~”. (4.28) 


Thus as0 > y > — 2, and X, < 0, we see that the Norm must decrease to zero. In the 
linear case (y = 0) this decrease is exponential. This proof is subject only to the proviso 
that At be small enough so that (4.1) will give us an approximation of the time derivative 
of the Norm that reproduces the sign properly. As the norm tends to zero, we see that 
our scheme is unconditionally stable, at least for this special class of problems. It should 
be noted that both the diffusion equation [2] and the radiation transport equation with 
power-law, temperature-dependent opacity belong to this class. 

Analogous arguments for the non-linear case corresponding to the arguments in- 
volving (4.10) to (4.12) indicate that the requirements for obtaining good time-dependent 
and steady-state solutions are about the same as the requirements for obtaining such 
solutions in the linear case. Therefore, the non-linear case has the same properties so far 
as stability and accuracy are concerned as the linear case. The only way in which the 
numerical solution to the non-linear case differs essentially from the linear case in rec- 
tangular coordinates is that an iterative procedure must be used on each time-step in 
the non-linear case. 

5. Extension to non-rectangular regions and three-dimensional problems with an 
axis of symmetry. Consider a simply-connected, two-dimensional region with more 
than one boundary point. Riemann’s theorem [5] implies that we may map it by a 
one-to-one conformal transformation into a rectangular region. This transformation 
amounts to a change of independent variables. If 


u= u(z, y) v = v(z, y) (5.1) 


is such a conformal mapping, then (3.1) is transformed [6] into 


97 a7 \ a) 
7(2 v + =) = h(u,v, wv) 1 , (5.2) 
Ou ol ot 


where 


—— au\? Ow : 

9 = (2) ° (2) 

(2): “ (*): (5.3) 
oY OY 


= g (u, Vv). 


II 
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Thus, defining 


H(u,v, WY) = h(u, v, ~)/g°(u, v). (5.4) 
we have 
== + * = h(x, y, W) o 
Ox Oy ot — 
(5.5) 
ay ay — Oy 
‘ ou" + ov" a Hu, Y, ¥) ot 


We may therefore first transform our region into a rectangle and then solve by the 
method of Sec. 3. It should be noted that when a non-rectangular region is to be solved, 
the iterative scheme of Sec. 3 is necessary, even in the linear case. 

One simple example of this method is a transformation to polar coordinates. Let 


u+ i = log (a + iy) (5.6a) 
or 

u = logr v= ¢@, (5.6b) 
where r and ¢ are the standard polar coordinate radius and angle. Equation (2.1) becomes 


0°06 ae ou «9 OF —" 
—“- —<. = € @ =e (5.7) 

Ou ov ot 
The singularity near the origin (u— — © asr—0) should be noted. This transformation 
actually carries a slit annulus into a rectangle. There are other more complicated trans- 
formations which will carry a whole circle into a rectangle. For instance, 


art+iu 


utw= | (1 — w)'” dw. (5.8) 


lor problems with an axis of symmetry, we adopt cylindrical coordinates as the basic 
point of departure. The Laplacian is 


3 oy + 1 i 4+ ay. (5.9) 


ror |r de dz” 


Viv = & + 


or” 


As we assume an axis of symmetry we may choose our coordinate system so y is in- 


dependent of ¢. Let us set 


a= ry, (5.10) 
Thus Eq. (3.1) becomes 
dw , Ww dw ow a 
or” + 027 nr, 2, «) at 4r’ (5.11) 


From the point of view of numerical solution we are free to add a term involving 
only w and independent of any derivatives, and the method of solution proposed in Sec. 
3 is still applicable. It should be noted that the iterative scheme will have to be modified 
to include this term, and if the coefficient of w at the advanced time contributed by the 
right-hand side of the difference equation analogue of (5.11) may change sign, we may 
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no longer pick 8 so that the error will change sign everywhere between the guessed and 
calculated values of a, but we must in some regions extrapolate for a new guess of w, 


9 


rather than interpolate throughout as we did in See. 3. 
The boundary condition to be applied along the axis of symmetry is that 


d(log w) , 5.12) 
ee ae co Pe 4 


d(log r) 
This condition ensures that y will be finite on the axis. 

For many problems with an axis of symmetry, cylindrical coordinates are not the 
most convenient. Other axial shapes may be treated by conformally mapping them 
into a rectangle. A sphere is conveniently treated by transforming (5.11) by (5.6) where 
we, of course, consider only the right half r — z plane (— 7/2 < v < 7/2). Thus for 
spherical coordinates we get 


I 


0 “wy 2 Ow Ww = 
a . 5.13) 


—- =e h(u,v,o) — — 
Ov dl | 


Ow 
au 
6. Numerical examples. In order to illustrate the method described above we 
programmed it for an IBM 704. Two sample calculations are described. 
We set up initial conditions corresponding to the solution of the linear problem, 


A(x, y, t) = sin (rx/I) sin (ry/l) exp (—2a‘r't/l’). 6.1) 


] 


We chose a’ At/(Ax)” = 1.02392228 so that the amplitude would be diminished by a 
factor of 10~*”’° at each time step. This choice enables us to check the accuracy at many 
times without excessive labor in calculating the analytic values. We ran this calculation 
with mesh points 15° apart (11 X 11 interior points), and advanced it 320 time steps. 
It was found that over the range of 20 decades through which the solution passed, the 
only discernible numerical error was a truncation error (the 704 does not round, but 
truncates instead) that caused the solution to decrease by an additional factor of 
(1 — 5.44 X 107°) at each time step. When a fixed number of decimal places are carried, 
the truncation error is expected to be proportional to a*(At)’/[(Ax)*(Ay)’]. The difference 
between the solution of the difference equation and the differential equation is that 
time flows 2.5 per cent more rapidly in the latter case. It should be noted that the explicit 
scheme is unstable for this case. The calculating time for this case was about 6 minutes, 
or about 10 milliseconds per cycle point. 

For our other example we chose radiation flow in a material medium. The equation is 


V'*(0') = 16K . 6.2) 
Ot } 
If we let y = K~*’°6*, (6.2) becomes 
' a 
Vy = 4y e. 6.3) 


[If there were a power-law, temperature-dependent opacity, the equation would be 


V-(6' 7 6) = 16K (6.4) 
Cc 
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By use of the identity 





e”" 7 oe =—— Vo" (6.5) 


n-+m 
and the transformation 
y = 6! (6.6) 


we can put (6.4) in the form of (3.1).] 

For our sample calculation we picked 6(Az)? = (Ay)? = — 1.5 and used 11 X 11 
interior mesh points. The other constants were chosen so that this is an optimal iteration, 
as defined by (3.3). For initial conditions we chose 


W(x, y) = 1 — sin (2r2/]) sin (2ry/1) (6.7) 


and we maintained the boundaries at unit temperature. We added 0.01 to y at the 
central point to prevent it from vanishing. Initially some 10 or 12 iterations were required 
per time step, but after 15 time steps, only about 2 iterations per time step were required 
for three-place accuracy. We make the initial guess of y at each time step by a linear 
extrapolation from the two previous times. The running time for this calculation was 
about 15 minutes for 120 time steps, or about 20 milliseconds per cycle point. From 
runs with other values of 6(Ax)° we estimate that the accuracy of this solution is about 
2 to 3 per cent. 

It should be noted that some care must be exercised to prevent the guessing of a 
negative temperature, as the occurrence of temperatures of different signs produces an 
instability which causes the solution to diverge. 
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BOOK REVIEWS 
(Continued from p. 360) 


Mathematik und Mechanik which is engaged, with the Association for Computing Machinery in the 
United States, in constructing an international algebraic language, to serve as a digital computer source 
language for mathematical problems. An agreement on terminology is a prerequiste for such an under- 
taking, and this work provides an acceptable basis, for translating programming terms between the 
English, French, Swedish, Dutch, and Russian languages. A cursory check revealed only a few omissions, 


e.g. object language, source language, mask, table look-up, assembly program. 


WALTER FREIBERGER 


Computational methods of linear algebra. By V. N. Faddeeva. Dover Publications, Inc., 
New York, 1959. x + 252 pp. $1.95. 


This textbook on numerical methods consists of three chapters: 
I. Basic material from linear algebra (62 pages 
II. Systems of linear equations (84 pages). 
III. The proper numbers and proper vectors of a matrix (96 pages). 
It is a translation of a noted Russian work and is unique in its clarity, elegance and scope. There is a 
wealth of useful numerical examples which, with the numerical tables, have been checked by the trans- 
lator. The emphasis is on methods useful with desk calculators. 

The first chapter presents the mathematical background, assuming knowledge of elementary deter- 
minant theory but not of matrices. Most proofs are shown and the author works in the complex field. 
The Jordan canonical form is given without proof, as is the theory of elementary divisors. There is a 
section on vector and matrix norms and on limits. 

The second chapter starts with Gauss’ elimination method in its various forms, recommends the 
square root method for symmetric systems and demonstrates partitioning and escalator methods. A 
welcome feature is a careful discussion of convergence criteria for the ordinary and the single-step 
(Seidel) iteration methods, with numerical illustrations, but gradient methods are omitted. 

The third chapter treats, firstly, of various ways of determining the characteristic polynomial, 
methods associated with the names of A. N. Krylov, Samuelson, Danilevskii, Leverrier, Faddeev, and 
of the escalator and interpol ition methods for obtaining the eigenvalues. The classical iteration method 
for finding the largest eigenvalue, with and without accelerated convergence, is described, with - 
differencing chosen as the method for determining the following eigenvalues. 

The author has unfortunately not included a discussion of how automatic computing equipment 
affects the choice of techniques. If she had, she would no doubt have presented Jacobi’s method for 
eigenvalue analysis. There is a bibliography of 40 titles and a useful index. The translation is readable 
and the book is invaluable for everyone with any occasion to analyse linear systems. 


WALTER FREIBERGER 


Operations research—methods and problems. By M. Sasieni, A. Yaspan, and L. Iriedman- 
John Wiley & Sons, Inc., New York, and Chapman & Hall, Ltd., London, 1959. 
xi + 316 pp. $10.25. 


The purpose of this book is largely to provide a collection of illustrative problems for use in courses 
in operations research. The problems (many of which are in the form of worked examples) are primarily 
applications of elementary probability theory; they are classified by background topic under chapter 
headings of sampling, inventory, replacement, waiting lines, competitive strategies, allocation, and 
sequencing. A short general discussion is given at the beginning of each chapter. 

The keynote of the book is simplicity, and it is refreshing to note the author’s statement that ‘‘Al- 


though a great deal of elegant mathematics appears in the literature, far more problems have been 


(Continued on p. 408) 
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AIRFOIL IN A SONIC SHEAR FLOW JET: 
A MIXED BOUNDARY VALUE PROBLEM FOR THE 
GENERALIZED TRICOMI EQUATION* 


BY 
C. C. CHANG anv T. S. LUNDGREN 


University of Minnesota 


Abstract. In this paper, small perturbations of a non-uniform two-dimensional 
flow of a compressible inviscid fluid are considered. It is shown that for a particular class 
of main stream Mach number distributions, which are characterized by a sonic line along 
the x-axis, the linearized shear flow equation may be transformed into the generalized 
Tricomi equation. The mixed boundary value problem which results from considering 
perturbations generated by a two-dimensional camber surface is formulated and solved 
by utilizing the Wiener-Hopf technique. 

1. Introduction. In recent years, Lighthill [1] Chang and Werner [2] and others 
have treated two-dimensional linearized compressible shear flow problems. In this 
paper a parallel but non-uniform stream of compressible fluid is perturbed by a thin 
two-dimensional airfoil which is located along a sonic line. The main stream Mach 
number distribution, 1/(y), resembling that approached by a sharp edged sonic jet 
under the influence of viscosity, is assumed symmetrical about the sonic line at y = 0. 

Let p(x, y) be the perturbation pressure, P be the constant freestream pressure, and 
a(x, y) be the slope of the streamlines in the disturbed flow. Also let M, = M(0). a and 
p, if small enough, can be determined from a function ¢(x, y) satisfying the linear equ- 
tion [3] 
a6 


Ox” 


x) 2 dM a¢ (1.1) 


(1 — M* 5 = 
, 7 oy” M dy dy 


oS 


¢ is related to p(x, y) and a(z, y) by 


————E— en _ 
ax (/2)PM2 = ©? (1.2) 
d@_ _. M’ 
oy — 2 M: a. (i .3) 


Here C,, is the local pressure coefficient and y the ratio of specific heats. It is generally 
accepted that (1.1) is a valid linearization of the equations of fluid motion provided that 
dM /dy is large enough at points where M = 1. This condition is not satisfied for the 
type of Mach number profile considered in this paper. It appears reasonable, however, 
that (1.1) is valid if higher derivatives of 17 are large enough at sonic points. 

Let the airfoil have unit chord length and be a camber surface so that the slopes 
of the upper and lower surfaces both equal a,(x). The condition that the flow be tangent 
to the airfoil can be obtained from (1.3) as 


Idb(x, O 
$9.0) a nhigit. ee 2s. (1.4) 


*Received September 3, 1958. 
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For subsonic flows, the disturbance must die out at infinity. This condition can be 
accomplished by imposing that ¢ approaches zero at infinity. The problem can be some- 
what simplified by specifying boundary data along the whole z-axis instead of only 
on the airfoil. It is shown in Appendix I that, as in incompressible potential flow ¢ has 
a jump discontinuity which is equal to the lift coefficient, C, . This jump can be taken 
conveniently to be along the positive z-axis behind the airfoil. In addition, because of 
the symmetry of the Mach number profile and the anti-symmetry of the airfoil, the 
transformation y — — y, ¢ — — ¢ leaves Eqs. (1.1) and (1.4) unchanged. This implies 
that ¢ is an odd function of y. This fact and the jump condition require 
¢(z, 0.) = 0 x<0O —- 
(1.5) 
= —C,,/2 >i. 
Equations (1.4) and (1.5) give mixed boundary data along the entire z-axis. 

In general, solutions of (1.1) are difficult to obtain with any arbitrarily given M(y). 
However, as will be shown in Part II, for a particular choice of M(y), Eq. (1.1) can be 
transformed to the form 

Prt. 2ece ree, (1.6) 
Ox Oo} 


where 7 is a non-negative constant and Y is a new independent variable 


Y= bf M’*(y) dy b = constant. 


) 
Eq. (1.6) is commonly known as the generalized Tricomi equation. This equation has 
been widely studied as shown in the bibliography to [4]. However, no solution is known 
for the above mixed boundary conditions. 

The present treatment will be divided into two parts. In Part I, the mathematical 
solution to this problem will be derived. In Part II, the application to an airfoil in a 
sonic shear flow jet and the details of the corresponding Mach number profile will be 
considered. 

Part I. SOLUTION OF THE GENERALIZED TRIcoMI EQUATION WITH 
Mrxep BounpDARrY CONDITIONS 

2. Statement of problems and procedure for solution. ‘The problem is to solve the 

generalized Tricomi equation 


- a "7 : O, Yr = 0. n=0 (2.1) 
Ox c 


with the mixed boundary conditions 


¢(z, 0) = 0, s <0; 
( C 
d(x, 0) = —5 . SS 3: (2.2) 
4] ( 
2, % = Oz), <2 < f: 
oY 


o(z, ~) = 0. 
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Here C is a constant and @(x) is bounded and integrable. If ¢(z, 0) were given along 
the whole z-axis, the problem would be a Dirichlet problem which could be solved easily 
by use of classical techniques such as the Fourier integral or two-sided Laplace trans- 
forms. By applying Wiener-Hopf techniques [5, 6] to the present problem, the missing 
Dirichlet data, ¢(2, 0) on the interval (0, 1), can be found. This is carried out in two 
steps. In See. 3 a related auxiliary problem is first solved, namely, the generalized 
Tricomi equation with d¢(x, 0)/dY specified on the positive z-axis, and ¢(z, 0) = 0 
on the negative z-axis’. This is done with the aid of two-sided Laplace transforms and 
the Wiener-Hopf technique. In the main problem 0¢(z, 0)/dY is given on (0, 1) and 
(x, 0) is specified on (1, ©), consequently, the auxiliary problem gives an integral 
equation for the unknown quantities ¢(x, 0) on (0, 1) and d¢(z, 0)/dY on (1,). In 
Sec. 4 this integral equation is solved by using Mellin transforms with the Wiener-Hopf 
technique. 

3. Solution for ¢(x, 0) on (0, ~ ) when d¢(zx, 0)/d Y is specified on (0, ~ ) and g(x, 0) = 


Q0on(—-,0). If¢(z, Y) is a function which satisfies 
=! a2 . 
ale nh FRO 228 (3.1) 
Ox oy 
and the boundary conditions 
d(x, 0) = 0, 2 <= @: 
Ad(a, 0) ' 
CAT, — O(2), 2 >O0; (3.2) 
0) 
g(r, ©) = 0, 
it will be shown that on the positive x-axis, the solution is 
(0,0) =r fe -9rmrae | adn - O°" Pda, 8.3) 
nt 4) Jo Jt 


where 
a —(1/n + 2)°"”"'*P(n + 1/m + 2)/TU/m + 2). 
First, let the boundary conditions (3.2) be restated as 


d(x, 0) = fi(x)H (2), 


— = 6,(z)H(2) + 6,(2)H(—2), (3.4) 


where /,(2) and @,(x) are unknown functions and H(z) is the Heaviside step function, i.e., 
H(z) = I, z>o0 
= (), $< @. 


1This problem differs from that of Carrier [7] in that Carrier’s boundary data are given on a line 


perpendicular to the parabolic (sonic) line. 
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Secondly, let 


nx 


G(s, Y) = | | e “d(x, Y) dx (3.5) 


be the two-sided Laplace transform of ¢(z, Y). Taking the Laplace transform of (1.1) 
and (3.4) gives 


re) “® 


= = 0 
sy? 


Ys’ + 
Cc 


0P(s, 0) 


—- @,(s). + 0,(s)_ , (3.6) 
oy r 


&(s, 0) = F,(s) 
Gig, ao) = |, 


where 0,(s), , 62(s)_ and F,(s), are defined by 


@,(s). 2 | e **6,(x) dx, 
0,(s)_ = | e *6(2) az 


F (3), = | etka) ae. 


Some imposed properties of these functions and their consequences are given below: 
(i) The given function @,(z) is assumed to approach zero like exp (— ev) asx — ©; 


thus 9,(s), will be analytic for Re(s) > — e where « > 0. Since 6,(x) is a bounded 
function, 9,(s), = 0|s|ass— o. 

(ii) 8,(s). , an unknown function which is to be determined, must approach zero 
algebraically as x — — © so that 6,(s)_ will exist for Re(s) = 0 and be analytic for 
Re(s) < 0. The assumption 6.(z) = O(— x)’, p < 1 as x — O assures that @,(x) Is 
integrable near the origin. This condition implies that 0,(s). = 0| s |~'*? as | s > ©, 

(iii) F,(s). is also an unknown function to be evaluated. The assumption here is 
that f,;(z) = 0 exp (— ex) asx — © and is integrable near the origin. These conditions 


1 


imply that F,(s). is analytic for Re(s) > — eand0|s|"*,¢g>Oas|s|— ©. 
The relevant solution of Eq. (3.6) is 


. 4 


where s must be purely imaginary to satisfy the condition 6— 0 as Y — o. Differentia- 
tion of (3.7) yields 


d® \ 1/2 un+l/2 1/8 2 3/2 
— = —A(—s’)'”Y"*'’K,, wil —s°)'* ——— yen (3.8 

dy ) s+iznra] (8) n+2 ) 
Evaluating (3.7) and (3.8) at Y = 0 and using the boundary conditions from (3.6) we 
can write 
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(1/n + 2) °°" *"(—s)'7"*" 
&(s, 0) = r/2 — - —innns A as f(s), , 
sin x(n + 1/n + 2), n + 2) (3.9) 
. \ +1/n+2) —1/2n+4 
d®(s, 0) _ = (1/n + _—™ eee (=8) A= 60. + Oe.. 








- ——T, 
d) atm (1/n + 2) 
Eliminating A from these two equations, we find, 


9, { + + 0,(s)- (3.10) 


Oe ia +2 ’ 


aF (s), = 


where a is given in Eq. (3.3). 

As it stands, Eq. (3.10) is valid only for purely imaginary s. In order to use the 
Wiener-Hepf technique, an expression which is analytic in a strip parallel to the imagi- 
nary axis is desired. Consider the function (— s*)'’"**, Re(s) = 0 which is real and positive. 
How should this function be continued analytically into the complex s-plane? Consider 
the product (— s)'’"**(s)'’"** for complex s. (— s)'“"** is defined to be real when s is real 
and negative and a cut is introduced in the s plane along the positive real axis so that 
this function will be single-valued. s'’"** is defined to be real when s is real and positive 
and a cut is introduced along the negative real axis. It is easily shown that when s is 
imaginary this product is equal to (— s°)'“"**. The above product is taken as the con- 


tinuation of (— *. New let a small separation of the branch points be made by 





‘.* D ‘ l/n+2 ‘ ne rr ° ° ° ye 
writing (— s (s + 6) , € > 0. See Fig. 1. This is the same e introduced in (i). 
Ss plane 
ain = 
cut tor (3+6)"3 — for (-s)"** 
€ 








Fig. 1. Analytic region of (—s)!/"*2(s + «)!/*?, 


If this function is used in place of (— s°)'’"*? in (3.10) we get 


aF (s).(s + )'"? = — oO ee (3.11) 


It is assumed that the solutions of this equation will satisfy (3.10) when ¢ goes to zero. 
This idea was used in [9]. Now the Weiner-Hopf technique can be used since both sides 
of Eq. (3.11) are analytic in the strip — « < Re(s) < 0. 

We want to express 9,(s),/(— s)-“"** as the sum of two functions, one analytic in 
the left half plane, the other analytic in the overlapping right half plane. By Cauchy’s 
Integral Formula [11], 


0, (s), 1 re 1 O,(w). 
- —s ] re 
ane - = ‘ Sar iy eo wa | 
: M1 Jy ,-io W s(- w) — 
1p _ 1 _@@), 
vi 2Qri i or s(— (—«)'” 3 dw, 


where — « < y, < Re(s) < y. < 0. The first integral in Eq. (3.12) is analytic 
<2, While the second is analytic for Re(s) > y, . Equation (3.12) can be 


for Re(s) Y2 
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written as 


sn = o,(s)_ — o,(s), , (3.13) 





where o2(s)_ , o:(s), are defined by Eq. (3.12). Using Eq. (3.13), Eq. (3.11) becomes 


; sips a 8,(s)_ : : . 

af,(s).(s + 6)'”"*? + o,(s), = 7_,)izare + 92(8)- = E(s). (3.14) 
—8) 

Since the two sides of this equation are analytic in overlapping half planes, they are 

sufficient to define an entire function E(s). That is to say, the right hand side is the 

analytic continuation of the left hand side. From the estimates made in (i), (ii), (iii) 

and since y, and y, are of the order of | s |~*, it is seen that 


E(s), = O(| 8 |°""?-"*9, 
E()- = O(| s |r"), 


as |s|— ©, where p, g < 1. The derivative of E is bounded and goes to zero for large s. 
Therefore, by Liouvilles’ theorem [11] F is a constant. But the above shows E(s)_ — 0 
as|s|— o. This implies that E = 0. This fact allows the evaluation of F,(s), . From 
Eq. (3.14), 

= o,(8) 


(s rs e): nmt+2 9 


Let « — 0 and substitute the integral expression for o,(s), . This gives 





Fi(s), = —a Re(s) > —e. (3.15) 


- 7 6,(w). 





oo ee ae | ie ida 
F,(s) = a's" Oni r= @ (=a)? dw, Re (s) > 0, (3.16) 
where 7; = O since e > y,; > 0. 
Taking the inverse transform of (3.16) we obtain 
~ ne” peer we O,(w). ” 
¢(z, 0) = f, (x) =a 'L {s sania as [ —— ) + as} (3.17) 
2m J_io 8 — w (—w) 


The inverse transform’ on the right hand side is the convolution of 


1/n+2)-—1 \ 
wnn £m) 
oe) Rie 
with 
u4 .f +3: as} « 80 = / gt? Ns dy me Biz 4la), 
Nt JI atin (—w) 


272 J-i0 8 — w (—w) 








where ¥(z) is defined by this equation. Use of the convolution integral (the inverse 
transform of the product of two transforms) gives 


f(z) = a” mcs f (2 — )"°/"*? “"Ye) dt, «= at > Oz (3.18) 


2All Laplace and Mellin transforms used in this paper can be found in [10]. 











1960] AIRFOIL IN A SONIC SHEAR FLOW JET 381 


But ¥(é) is the convolution of 


joo 


oj | eile), do = 0,@H®) 





with 
- aie ~ aie a! (—#)'/"*?H(—2z) 
Qmi | “a du = T(i/n + 2) 
Therefore, 
/e\ ] ‘ \ (1/n+2)—1 
L(t — —— _— 
vie) Tare i, Ke é) dn, §>0. (3.19) 


Substitution of Eq. (3.19) into Eq. (3.18) yields Eq. (3.3) as was to be shown. 
4. The case where 0¢(z, 0)/dy is given on (0, 1). The result of Sec. 3 is 


a” [ (1/n+2)-1 = (1/n+2)-1 
(r+. 0) = a ey 2 _ atlas 
g(x, rane |, &-9 at | A:(n)(n — &) dn, (4.1) 
o<s< @, 
where 6,(x) = d(x, 0)/dy is supposed to be known on 0 < x < o. In the problem 


stated in Sec. 2, d(x, 0)/dy = 6(x) is specified only for the interval (0, 1) while ¢(z, 0) 
is specified for x > 1. Hence Eq. (4.1) is an integral equation for the unknown functions. 
It will be shown that if 


als, 0) = 6(2) = 0(2)H(1 — z) + 0,(2)H(z—-1), 2>0 (4.2) 


and 


d(x, 0) = f(x)H(1 — x) - CH — 1), x>0 


where f(x) and @,(x) are unknown functions, then the desired solution of Eq. (4.1) is 





C r@s/a.+:2) ff” ccrtjeces= ae 
é(7, 0) = f(z) = 3 PC * +9) I a ) 1] ae ind ) dt . . 
"ee taeda r (1/n42)-1 . (1/n+2)-1_1 ) 
= Sa ad n+2)— ¢72 n+2 == n+2)- /n+2 
+p | €- aren de [| — Verna 0(a) da, 


when z is between 0 and 1. 

To prove Eq. (4.3) the same procedure is followed as in Sec. 3, using the Mellin 
transform in place of the Laplace transform. Operating on both sides of (4.1) with the 
Mellin transform, defined by 


M,{f@)] = [ 2*""f(a) dx, 
yields the relation® 


rfl — s — (1/n + 2)]T[s + (1/n + 2)] pr 
rd — §)P[s + (2/n + 2)) M,[x*”""*6,(2)]. (4.4) 





aM, |¢(x, 0)] = 


*This is the Mellin transform of a repeated fractional integral. See [10] , vol. 2, p. 183. 
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A similar operation on Eq. (4.2) yields 
, ' ' : C ?— C I 
M,|¢(x, 0)] = M,[f(x)H — x)] - 5 M.[H( — 1)) = F®. + rege 
“ zs 
M, [x*’"*?6,(x)] = M,[2?’"*?0(2)H(1 — x)] + M,[2?’"*? 6,(2)H (a — 1)] (4.5) 
= T(s). + T;3(s)-_ , 
where 
al 
F(s). = | x "f(x) dx, 
] 
— = —M,/|H(z — 1)], Re (s) < 0, 
al 
ie, = | 2s Ox) dz, 
T.(s)_ = | ea G(x) dx. 
Substitution of Eq. (4.5) into Eq. (4.4) yields 
; C1 _ Ml—s—(imt+2)Tetidmt2)) on. . mse 
ars). + ¢ 2 3. 7 (il — s) rs + (2/n + 2)] Ne)» + Tals). i 


As in Sec. 3, the regions of analyticity and the asymptotic behavior of F, 7, 7; are 


required and are shown as follows: 

(i) F(s). is an unknown function. It is assumed f(z) = 0(2'“"**) as  — 0; then F(s), 
is analytic for Re(s) > — 1/n + 2. In order to determine the asymptotic behavior of 
F(s). , substitute a exp (— 2) into the defining integral above, so that it takes the 
[7 p<lasxr—1 


form of a Laplace integral. It is easily seen that if f(z) = 0|1—a 

then F(s), = 0| s {~*~ as| s|— o. This assumption is necessary in order that F(s), 

exist. 
(ii) T(s), is a known function. Since 6(x) is a known bounded function, 7'(s), 


analytic for Re(s) > — 2/n + 2. The substitution x = exp (— z) shows that 7(s) = 


Olsltas!s|— oa. 

(iii) 7(s)_ is an unknown function. By comparison with known potential solutions 
(i.e., potential flow about a flat plate with circulation) it is assumed that the unknown 
“"*") asx — © so that 7,(s)_ is analytic for Re(s) < 0. If it is 


is 


function 6,(z) = O(x 
assumed that @;(z) = O(a — 1)"*,g < lasx-—~1, the substitution z = e’ obtains 7;(s)_ = 
O|s|"'**as|s —> 

(iv) The asymptot 
found by substituting : 


ic behavior of ratios of gamma functions in Eq. (4.6) may be 


exp (— ¢) into the beta function representation, 


I'(s) | [ s- 1] .? 1 Bae 
—— 22 6 — J ( 
r(s + ») Ti) Jp a r) oe 


From this it is found that the function I'[s + (2/n + 2)]/T[s + (1/n + 2)] which is 


analytic for Re(s) < 1 — (1/n + 2) is of the order | 8 |’ a s —> ©, 
while I'[1 — s — (1/n + 2)]/T(1 — 8) is analytic for Re(s) < 1 — (1/n + 2) and is of 
the order S | ee as S — Oo, 
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Rearrangement of Eq. (4.6) gives 
r[s + (2/n + 2)] Tl —s — (1fn + 2)] 











al’(s), Tis (Un r 9 i Ta —=»s) T';(s)- 47 
_ C1il®rfs+Qnr+2)) ,, rl —s — (1m+2)] 
1>5.Tl+dnm+2)) 7 2 rd — s) 


A little reflection on (i), (ii), (iii) and (iv) shows that all the functions occuring in Eq. 





(4.7) are analytic in the strip — (1/n + 2) < Re(s) < 0. As in Sec. 3, we can write 
fl —s — (I = — —w- 
ri =~s—-(m+%B my _ 1 [ _1_ Til-»o-(m+2))] x, sl 
(1 — s) 271 Jg.-i10 Ww — § (il — w) 
J oan l —w-— (] oa 
oer f Se Se ee ee 
271 Js.-io W—S rd — w) 
where — (1/n + 2) < 8, < Re(s) < 8, < 0 and \,(s)_ and X,(s), represent the first 


and the second integral on the right hand side. \.(s)_ is analytic for Re(s) < 6, while 
\, (s). is analytic for Re(s) > 6, . Also, by subtracting the principal part of the singu- 
larity at s = 0 

1 Ti[s + (2/n + 2)] _ JI Tis + (2/n + 2)] _17T in + Bh 1 T2/n + 2) (4.9) 
s T{fs+ ( ] ls T'[s + (1/n + 2)] T(1/n + 2) s. TU/n + 2) en 
The bracketed term on the right is analytic for Re(s) > — 2/n + 2, while the second 
term on the right is analytic for Re(s) < 0. Introducing the results of Eq. (4.8) and 
Eq. (4.9) into Eq. (4.7) gives 

ris + (2 n + 2)] 


l 
ie s 








A) = a9)* Tie + (1/n + 2)] + 
C {1 Tis + (2/n+2)]_ 1TQn+2 
+4515 Tis (Bu 2] _ 1 Tien + Bi) aa 
ri —s — (Ilm+2)] _ £TQ@n+2)1 
r(l — s) T(s)- + dals)- 79 TU /n + 2) s_ 


As in See. 3, since the left and right sides of this equation are analytic in overlapping 
half planes, an entire function £,(s) can be defined. From the estimates in (i), (ii), (iii) 
and (iv) it is found that 

E (s),. — 0 | 8 {(1/n+2)—-l+p 

” os (4.11) 

and 

E,(s)_ = _ O|s |~ (1/n+2)—-l+@¢ 
as 

|s|—o, where p,q <1. 

These estimates imply £,(s) = 0. This determines F(s), and 7'(s)_ in Eq. (4.10) as 
follows 


‘[s + (1/n + 2)] 


1 





s + (2/n + 2)] di(8)+ (4.12) 


T 
T{: 
fa it (2/n + 2) Ps + (1 in + 2H) Re (s) > —1/n+ 2 
~ 2 |g s T(1/n + 2) I'[s + (2/n + 2)] ? 


F(s), = —a” 
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il — s) 


eee | St ee 
Pils rl —s— (im+2)]” - 
1.13) 
C T(2/n + 2) 1 r(l —s Re ( - () 
t= ——$—— »(s) < 
12 Tin +2)sTil—s—(Unr+2] “ 


The next step is to take the inverse Mellin transform of Eq. (4.12) using the Mellin 
convolution in the same way that we used the Laplace convolution in Sec. 3. The result 
is Eq. (4.3) which was to be demonstrated. 

5. Solution of the mixed boundary value problem stated in Sec. 1. The original 
problem to be solved was a mixed boundary value problem. Since from Sec. 4, d(x, 0) 
is known on the interval (0, 1), we know ¢(z, 0) for the infinite interval. This reduces 
the mixed boundary value problem to a much simpler Dirichlet problem. 


Od . Oo 


y — 0, Y 2 0, n = 0; 
Ox oO} E ee 
».1) 
Ov 0 Do\2 = 6 ce m- < x 
C c } 0. 
Here 
0 r < 
P ? <= 2 I 
~ C/2 > 1 
where f(a) is defined by Eq. (4.3). Use of Laplace transforms as in Eq. (3.6) yields the 
bounded solution shown in Eq. (3.7). From Eq. (3.9) 
- = 1/2n+4 
¢ ¢ ( = - ca gas —s : A = pm S 5 2) 
2(1/n + 2) © sin (r/n + 2)Tin + 1/n + 2) 


where %(s) is the two-sided Laplace integral of ¢,(x). Solving Eq. (5.2) for A and 
substituting it into Eq. (3.7), 


9) : , 
®@ = -(l/n + 2 sin (r/n + 2)T(n + 1/n + 2)(—8 
: 5.3) 
' 2 
. - 7 ———— J , Re(s) = 0 
wT 2 
The inversion of this gives the result 
3 b(t) Y dé 
d(x, } d, — —— - 
J-w (te — &) + (2/n + 2)} 
where d, = 7°" (2/n + 2)° ssinlaa T{(1 nt+2)+] 2) T'(n + I1/n + 2) sin (7r/n + 2). 


Since this result will not be used in Part II, no further comment will be made. 


Part II. AtrRFoIL IN A Sonic SHEAR IF Low JET 


6. Solution of shear flow problem. As seen in Sec. 1, for a symmetrical Mach 
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number distribution the potential ¢ must satisfy 


d¢ , 96 _ 2 dM dg _ (6.1) 


(1 — M’) 


ax” * dy” M dy dy — 


for y = O, and the boundary conditions 


o(2, 0) = 0 x < 0, 
= —C,, 7 > i, (6.2) 
Te os eld ee <4, 
oY 
g(x, ©) = 0. 


Introducing the new independent variable 
fas / M*(y) dy, (6.3) 
0 


where } is a constant, Eq. (7.1) is transformed into 


1— M'd¢ , 06 _ 
b°M* dx” ” 0. (6.4) 


Let 1/(y), an arbitrary function until now, be defined by 


_ m| 2 P vy . 
ae = y"= (2 [ M? ay) (6.5) 


where 7 is a non-negative number, so that Eq. (6.4) becomes the generalized Tricomi 
equation treated in Part I. The solution of (6.5), which gives the Mach number dis- 
tribution, will be investigated in Sec. 7. The transformation by Eq. (6.3) also changes 
the third equation in the boundary conditions, Eq. (6.2) to 


. j ®., fae 

cee _ — ta) O<a<1. (6.6) 
It is seen immediately that this is the problem solved in Part I with the notational 
change 6(x) = — 26)(x)/bMs and C = C, . The solution of this problem involves the 
constant C, which has been presumed specified in advance. To determine C, the Kutta- 
Joukowsky condition should be used. This condition states that the flow must leave the 
trailing edge smoothly. In a linearized approach, this means that the streamline leaving 
the trailing edge of the airfoil should have finite slope at the trailing edge. In Appendix 
II it is shown that the slope is finite provided that 


2r(1/fn +2) (I1m+2)”""_ f' 


r(2/n + 2) Tin + 1/fn + 2) Jo 





C;, = (1 — 9) "*?-*9'/™*?[—Qa(n)/bM3] dn. (6.7) 


This gives the lift coefficient in terms of the airfoil slope. The pressure coefficient on the 
airfoil can be calculated by taking the z derivative of Eq. (4.3), since C,(z, 0) = 
d¢(x, 0)/ax. This may be accomplished more easily if the integration variable in the 


second integral in Eq. (4.3) is changed by the substitution z = 2/. The result is 
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; . : ; 
C,(z, 0) = e, | C—O a 
J1 Ou (6.8) 


at/2 - (1/n+2 l 
| (2 <a *) n "**1—2Qao(n) bM;] dn 


where 
a. 
™ ~ Tin + 1/fn + 2)T(/n + 2) 
It will be shown in Sec. 8 that Eq. (6.8) can be integrated with a simple result if ao(z) 
is a polynomial. 

7. Mach number profile. Certain properties of the upstream Mach number distribu- 
tion, M(y), which are defined by the solution of Eq. (6.5), can easily be seen without 
completely solving this equation. These will be enumerated below. 

(i) If nm = 0, M = M, is constant and less than one. In this case b(n = 0 by can 


2 


be expressed in terms of M, as bh = (1 — Mj)'”’/M; . Noticing how by occurs in Eq. 
(6.8), it is seen that b, is essentially the Prandtl-Glauert correction factor for incom- 
pressible flow. 

(ii) Forn = 0, M—1lasy—O0 while M —>Oasy— ~. 

(iii) To see how M behaves near y = 0, rewrite Eq. (6.5) as 


pm "MP dy = (2+*) (7.1) 


70 / 


and substitute a series, M = 1 — y’ (a, + ayy + ---), for MW. Solving the resulting 


equation for a, and » it is found that 


M =1-— }b"7y"+---. 7.2) 


This shows that in addition to being sonic at y = 0, the profile becomes increasingly 
flat near y = 0 as n increases. This is of some interest since, as is well known, purely 
sonic uniform flow cannot be linearized. Notice also that b — 0 implies uniform sonic 
flow, while from Eq. (6.8) the pressure becomes infinite as b — 0, so that the linearization 
has become invalid. 

(iv) The behavior for small n can also be seen from Eq. (7.2). If 0 < n < 1, M has 
a cusp at y = 0 which becomes sharper as n decreases. Upon taking the limit as n — 0 
the cusp collapses on itself leaving uniform subsonic flow. 

Since the solution of Eq. (7.1) for M as a function of y is difficult to obtain explicitly, 
the inverse of this function y = y(M) will be considered. Denoting M’ by é and differen- 
tiating Eq. (7.1) with respect to y we get 


l/n + 
pete d (154) dé 


dé 


Upon separating variables and integrating, 


n+2)/n i l d ] aw, t\' 
, “ee | E dé ( rs ) dg 


is obtained. An integration by parts gives 


an inA=-wy £3 254)" : 
b n « = —sz qua -_—. — -—> - =: f. (7.3 
Y= ( uM ) I. e\ * oo 
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This is the desired inverse function. If n equals 1 or 2 the integration can be carried out 
explicitly. For large n a series in 1/n is appropriate. The results of these calculations are 
given in Fig. 2. In this figure M is plotted versus yb“*””. If b = 1 the various curves 
are in the correct proportions as plotted, however, if b ¥ 1 the scale changes for each 
curve. If n is very large the profiles are again approximately in the proper proportions. 
Note that as n — © the profile approaches a uniform sonic jet of width 2/b; twice the 
reciprocal of 6 can be considered roughly to be an effective jet width. 














1.0 
M 
0.8+ nei 
n=2 
nid 
0.6 + ~ 
| fis 25 
! 
| 
0.4+ | 
| 
| n 2100 
' 
02+ | 
n #1000 
x3 
0 + + + + i + + 2 — 
(@) 0.2 0.4 0.6 0.8 1.0 1.2 4 1.6 


Fic. 2. Mach number versus vertical distance with n as parameter. 


8. Calculations of C,, C, and C.P. In this section the calculation of the pressure 
coefficient, lift coefficient and center of pressure (C.P.) will be carried out for the case 
where a (x) is a polynomial. In this case — 2ao(x)/bM% can be expressed as a sum of 
terms like a;x’ and the lift coefficient, Eq. (6.7) will be a sum of terms like 





T(1/n + 2)(1/n + 2)""*? [ (1/n+2)-1_ 5+ (1/n42) 
C.9* = : ' 
75 T(2/n + 2)I(n + 1/n + 2) Jo (I n) od dn 


This integral is a beta function which is easily expressed as 


TU /m + 2) fn + 2)""*? Td /m + 2)r[(ifn + 2) + jf + 1) (8.1) 
*T(2/n + 2)T(n + 1/n + 2) r[(2/n + 2) + 7+ 1] ° 


The pressure coefficient, Eq. (6.8), on the airfoil will be a sum of terms like 


p - oe a z/s x (1/n+2)=1 m 
Cori = Dn | fi — g) /a42)—1, (int ) as [ amy yi*! n+2) dn, 
Jz Jo 





C1; = 2a 


where 


—— _a(1/fn + ove , 
I = ~Tin + 1/n + 2)T0/n + 2) 


The inner integral, a fractional integra! of the Riemann-Liouville type can be evaluated 
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rij +1+(1/n + 2)] (ap /z)it 27m429 


pz/2 1/n+2)—1 
= — rare dn = T(1 2 = x/2 
i (: ") " ”" Un + 2) Te 71 + Qn + 2)] ** 


Then 
1 
Cy; = k, | —<_ so 
where 
h a(1/n + 2)”"*? T[j +1 + (1/n + 2)) 5-4, 
, = = y 2 
Tin + lfm + 2) Tlj + 1 + (2/n + 2)] 
By a change of the variable of integration, t = (1 — z)/z this becomes 
C I/n + 2)”"** Pij+1+ (1/n + 2)] 
= —/ r a - os ae r 
pi ’ : 9 nf. (9 a a7 
Tin + 1/n + 2) T[j + 1 + (2/n + 2)] (8.2) 


| gitar~tys 4. A! dit: 


The integration here can be easily carried out for any value of j by expanding (7 + 2)’. 
Let a moment coefficient be defined by 


“# =_- —9? | a( - dx. (8.3) 


This can be computed after C, is known. A more useful quantity which will be used is 
the center of pressure defined by 

CP. = C,/C, . (8.4) 

These quantities (C,; , C, , C. P.) are calculated below for the special cases of flat 


plate and circular arc airfoils. 
(i) Flat plate airfoil at angle of attack 6. 


In this case, using the above notation, a(x) = — 6 is a constant. Then 
— 2a,(x) 26 
7 > 2 = 2 = a 
bM; bM, 


From Eq. (8.1), 
I(1/n + 2) 


no , 28 (1 ) a t 2) 
’~ bMS \n + 2 r°(2/n + 2)0(n + 1/n + 2) 








(8.5) 


If nm = 0 (uniform subsonic flow) this gives the familiar result 


CG, = 2rh _ rh 
ai.) 6 6G wr” 

where b = (1 — M’)'”?/M”’, M constant. If n — @ (sonic jet) the result is C; = 86/b = 
4 hé. Here the jet thickness h = 2/b has been used. Equation (8.5) is plotted in Fig. ¢ 
From Eq. (8.2) the pressure coefficient 


- 


T(2/n + 2) er eee ren 
C, = oe ee ee (8.6) 


ee = -— y 
. (1 n+ 2)° 
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Fic. 3. Lift coefficient on flat plate versus n. 
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Fic. 4. Pressure coefficient on flat plate airfoil versus distance along airfoil. 
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Fic. 5. Center of pressure on flat plate airfoil versus n. 
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Fia. 6. Lift coefficient on circular arc airfoil versus n. 
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Fic. 7. Pressure coefficient on circular arc airfoil versus distance along airfoil. 
is obtained. In Fig. 4, — 2 C,/C; has been plotted versus x so that the “area” under 


each of the curves is unity. Notice that for n = 0 the leading edge singularity is like 
-1/2 . . . ° is 
xz **, As n becomes larger this singularity becumes stronger, approaching x 


The center of pressure is obtained by substituting Eq. (8.6) into Eq. (8.3). The result 


‘as n> ©. 


is the simple expression 
‘ae l _ 
Ca. oie, (8.7) 
This shows that for n = 0 the center of pressure is at the quarter chord and approaches 


the leading edge as a n increases. This is presented in lig. 5. 
(ii) Circular are airfoil at zero geometric angle of attack. Yor the circular are ao(x) = 


5(1 — 2x) where 6 is the angle which the leading edge makes with the free stream. ‘Then 
— 2a,(x) —26 46 


= + Qt, Ay = bM? ] = bM2- 


The lift coefficient and pressure coefficient are easily calculated from Eqs. (8.1) and 
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(8.2). We obtain 








26 - r(1/n + 2) 
C & 5) 9 n/n - [hae aS. 4 4 
ATE (n + 2/n + 4)(1/n + 2) rQn+2rn+iny+2)’ (8.8) 
o. a Le n + a .\l/n+2 
C-~ —-Pam +2 (n + 4)a'” x) , (8.9) 
The center of pressure is at x = .5 for all nm as can be seen from the symmetry of the 


pressure coefficient. Equations (8.8) and (8.9) are presented in Figs. 6 and 7. Note in 
Fig. 7 that asm — © the pressure coefficient becomes constant along the airfoil. 

9. Discussion. It has been shown that for a particular class of Mach number pro- 
files which are characterized by a sonic line, the linearized compressible shear flow 
equation can be transformed into the generalized Tricomi equation. The boundary 
value problem, which describes the perturbation of the main stream flow by an arbitrary 
two-dimensional camber surface, has been formulated and solved. The question of 
agreement with reality remains. 

It is probable that for large values of n the linearization breaks down unless re- 
strictions are made on the parameter b, which should be large as n becomes large. That 

, the effective jet width, 2/b, must be much less than the chord length. This conjecture 
is pets «1 by the following observation; by the foregoing theory the lift on an airfoil 
in the limiting ease m — © is the same as that calculated by a simple momentum ex- 
change brought about by the sonic jet turning through an angle equal to the trailing 
edge angle of the airfoil. Such a situation can only be realized if the jet width is much 
less than the chord length. 


APPENDIX I. 


Jump condition. Let C denote the curve indicated in Fig. 8, consisting of two 
straight lines and a circle of radius R. If V¢ has continuous derivatives inside C we can 


write 
[9 44 
| (28 dx + as dy) = 0 
J- \OX oY 
or 
oR a2 ao 
| (C,(@, 0.) — C2, 0] de + | 59 18 = 0. 
+ 
R c 
e 








Fic. 8. 
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Here 0¢/dx = C, has been used. C, must be continuous except across the airfoil itself, 


therefore C,(x, 0.) — C,(x, 0.) must be zero for x > 1. Then the first integral is the 
negative of the lift coefficient, 

: 
; ; , ; Lift 
C, = [ [C,(z, 0-) — C,(a, 0,)] dx = ss D)PM? 


/0 





while the second integral is ¢(R, 0_) — ¢(R, 0.). Hence 
C, = o(R, 0_) — ¢(R, 0,). 
Since C, is independent of R the jump across the z-axis must be constant. 
APPENDIX II. 


Kutta-Joukowsky condition. It is necessary that d¢(x, 0)/8Y be bounded as 
x — 1 from the right. That is, @,(a2) in Eq. (4.2) should be bounded. This condition will 
be satisfied only for a certain value of C. @;(x) can be found by inverting Eq. (4.13). 


Proceeding as with the inversion of Eq. (4.12) yields, after an integration by parts. 


8.(x) oT es a at’ r(2 _ 
Me) = SS : eee pena — — [(2/n + 2) 
Tin + 1/n + 2)T(1/n + 2) 2 


= | (l— 7 n “6(n) in| + eS oe 
J Tin + 1/n + 2)P(1 n+ 2) 
0z 1 
‘ | £ "(a — §) ° dé (§ — ny)” : °n iia n) dn. 


Therefore @,(z) will be bounded at x = 1 if the coefficient of (x — 1)~'”"** is zero, that 


is, if 
c a’ [ | /n+2)—1 1/n+2 a7 yg 
— ———— (1 — n) 1 (n) an. 
2 T(2/n + 2). . ’ =? 
Using the definition of a, Eq. (3.3), gives 
C i r(1/n + 2 e arene 1/t 
= (1/n + 2) ; a | (1 — )'’"*? "yn “A(n) dy. 
2 r(2/n + 2)0(n + 1/n + 2) sili 
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ON THE APPLICATION OF ELLIPTIC FUNCTIONS IN NON-LINEAR 
FORCED OSCILLATIONS* 


BY 
C. S. HSU 


University of California, Berkeley, California 


1. Introduction. A system having an equation of motion of the Duffing type is 
known to respond in free vibrations in the form of elliptic functions, and the period is 
expressible in terms of complete elliptic integrals. It seems appropriate, therefore, to 
analyze the response of such systems when they are subject to periodic forcing functions 
which are themselves elliptic functions. It is the purpose of this paper to examine this 
problem in detail. 

We shall show that there exist quite general cases where exact single-term and 
multiple-term elliptic responses can arise under forcing functions consisting either of a 
single elliptic function, or of a combination of these. It will be shown that all known 
results of the perturbation or iteration method are simply special cases of these exact 
solutions, and the former emerge easily from the present analysis if one expands the 
results of this paper in terms of the moduli of the elliptic functions and retains only 
terms which are linear in the second power of the modulus. Higher order approximations 
may easily be found from the multiple-term exact solutions by comparing orders of 
magnitude of the various components in an appropriate fashion. Finally, we shall 
demonstrate exact subharmonic solutions which are elliptic functions and which may, 
or may not, be simple. The simple elliptic subharmonic solutions reduce to the simple 
subharmonies defined by Rosenberg [6]** if one sets the modulus equal to zero. 

In view of the fact that there is a fair amount of uniformity concerning the notation 
used in connection with elliptic functions, we have adopted the commonly accepted 
notation of Whittaker and Watson [9]. 

2. Non-linear free vibrations. The equation of free vibrations of a Duffing system is 


m dz + k,.xz + k,,2* = 0. (1) 
dt 


Introducing r = ¢(k,,/m)'” and 6 = k,;/k,, , we obtain a new equation of motion 


dx 
7d? +-2+ Br* = 0. (2) 
aT 
Let us assume that the peak amplitude of the oscillation x is equal to A, and it occurs 
at 7 0. Then it is well known [9] that the solution of (2) is given by the Jacobian 
elliptic function en u. The exact form of the solution is 
x A en (pr, k), p = (1+ BA’)!, k = [gA?/2(1 + BA”)]}, (3) 


where k is the modulus of the elliptic function and p may be taken as the “circular 
frequency.” 

*Received September 16, 1958; revised manuscript received January 23, 1959. 

**Numbers in square brackets refer to the bibliography at the end of the paper. 
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If we have a soft non-linear spring, i.e., 6 is negative, the solution may be transformed 
[1] into the following form in order to have the modulus remaining in the usual range 
from 0 to 1, [2, 5], 
r , 2 soy\4 2s 2 } 
z= Asn [p,7 + K(k,), ky], p, = (1 + BA‘/2)’, k, = [-—BA°/(2 + BA’)]’. (4) 
Here K(k) is the complete elliptic integral of the first kind. 
Similarly, when 8A” < — 1, the solution can be put in the following form. 


z = Ans [p.r + K(k), ko], po = (—8A?/2)'”, = ke = [—(2 + BA”) /B AJ”. (5) 


We note that in this case the solution is no longer bounded. The meaning is quite 


obvious physically. When 8 is negative and A is so large that BA” < — 1, the total 
spring force—i.e., the linear and the non-linear part together—becomes negative. Under 
this circumstance the mass will move away from the origin x = 0 indefinitely. As a 


matter of fact, the quantity A should be interpreted here as the amplitude at 7 = 0, 
not as the peak amplitude. 

Both the functions enu and snu have a period in u equal to 4K. The complete elliptic 
integral of the first kind, K, is then the quarter period of the free oscillation. The period 


in 7 is given by 


T, = 4K/p. (6) 
For the case of a hard spring it is 
T, = 411 — 2k’)'’K(k). (7 
For a system with a soft spring it is 
i 11 +h K(k,). (8 


Figures 1 and 2 indicate these relationships between 8A°, k”, kj and T, in graphic form. 
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Fic. 1. Free vibration period and amplitude as functions of the modulus for a system with a hard spring. 
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Fig. 2. Free vibration period and amplitude as functions of the modulus for a system with a soft spring. 








When the non-linearity and the peak amplitude are such that | 8A” | is small in 
comparison with unity, the modulus k will also be small. Expanding the quarter period 
K in powers of k* and substituting the expansion into (7) we find for the period 


T, = 2n[1 — 3k’ + OK). (9) 
The equivalent circular frequency p,, can be easily found as, 
p., = (2r/T,)? = 1 + 368A’ + 0(8°A’'). (10) 


This agrees with the result obtained by the perturbation method when | @ | is assumed 
small [8]. 

3. Some exact solutions of forced vibration. In this section we shall examine the 
response of a Duffing system to forced vibration. We assume that small quantities of 
positive damping are present in the physical system, but absent from the equation of 
motion. Consequently the periodic solutions of interest are those usually referred to as 
steady state vibrations. We shall assume here that the effect of small damping on the 
motion is slight. 

Let us consider a system whose equation of motion is 


“ +2+ 62 = F(n), (11) 
where /'(r) is the external forcing function divided by k,, . The forcing function F(r) 
in (11) is usually taken to be a simple harmonic function of an arbitrary circular fre- 
quency. Here, however, this is not the case. Instead we shall try to find out what type 
of external forcing function will excite a steady state response, which is an exact solution 
of the equation of motion, and which is expressible in a simple form. 








396 C. S. HSU [Vol. XVII, No. 4 


3.1 Simple elliptic response. First let us ask under what condition the response 
x as a function of 7 will be proportional to the forcing function, leaving the precise form 
of the function undetermined. Putting P(r) = Bz, the equation of motion becomes 


a + (1 — B)x + Bx’ = 0. (12) 
The solution of this equation can be readily found to be 
x = A cn(pr,h), (13) 
where 
p= (1—B+ 8A’), k = [8A°/211 — B+ BA’)]'”. (13a) 


We note that there are four parameters, A, B, p, and k involved in the solution and 
two relationships (13a) between them. Hence, any two of them may be arbitrarily 
chosen and the other two will be determined by (13a). Such a set of parameters will 
assure a response like (13). For example, if we take p and k as given we have 
A = (2/8)'’*pk, B=1-—p(l — 2k). (14) 
The response x and the necessary forcing function F(r) are respectively 
xr = A cn(pr,k) = (2/8)'*pk en (pz, k), (15) 
F(r) = BA en (pr, k) = (2/8)'*pk{1 — p°(1 — 2k°)] en (pr, &). (16) 
The parameters p and k in en (pr, k) play important geometrical roles. The modulus k 
controls the proportions of the higher harmonics in the Fourier analysis while p and k 
together determine the periodicity of the function. The period in 7 is given by (6). 
Figure 3 shows curves of 7’, as functions of p and k. 
A response consisting of only a single term of elliptic function like (13) will be called 
a simple elliptic response. It should be emphasized here that not every elliptic forcing 
excites a simple elliptic response. The amplitude of the elliptic forcing function must 
bear a definite relationship to the parameters p and k. We shall denote those forcing 
functions satisfying this relationship as favored elliptic forcing functions. 
In case of a soft non-linear spring, 8 < 0, the solution of (12) is 


x = Asn (p,7, k,), 
where 
Dp 1 — B+ BA*/2)'”, k, [—8A°/2(1 — B + BA’/2)]'”. (17) 
If we consider p, and k, as given we must have 
B 1 — pi(1 + &)), A = (—2/8) “p k, (18) 
and the forcing function 
Fir BA sn (p,7, ki) = (—2/8)'"p,k, [1 — pi(l + k)] sn (p17, k,). (19) 


Equations (17) may also be obtained through the application of the first order 
transformation mentioned in Sec. 2. From here on we shall not, therefore, treat the cases 
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Fic. 3. Period in 7 of the elliptic functions cn(pr, k) and sn(pr, *). 


with soft springs separately. Whenever the modulus k appears outside the standard 
region of 0 to 1, one of the first order transformations must be used.* 

The quantity 1/B is analogous to the well-known magnification factor—i.e., the 
ratio of response amplitude to forcing amplitude. In the linear case it gives the response 
amplitude for a forcing amplitude of unity. Near the point where 1/B approaches 
infinity resonance occurs. With nonelinear systems it is better to take a somewhat 
different view. In linear cases it is customary to consider the forcing amplitude as given 
and the problem is to find the response amplitudes. In the non-linear case the response 
amplitude is intimately connected with the period, which is determined by p and k, 
and the wave form, which is determined by k. It is therefore desirable to consider the 
response amplitude as known and the problem is to find the necessary forcing amplitude. 
Hence, instead of 1/B we shall discuss the factor B, which is the ratio of forcing amplitude 
to response amplitude and will be called the forcing amplitude factor. Figures 4 and 5 
show the curves of B vs. p with k* as a parameter for hard and soft non-linear springs 
respectively. 

In Figs. 4 and 5 the curves of k* = 0 correspond to the usual linear case. B is given 
by (1 — p’) or the magnification factor 1/B = 1/(1 — p’). Each of the curves in Fig. 
*In a slightly different form, the simple elliptic response under one term elliptic forcing discussed 
here was treated by Helfenstein [3]. The material presented here was retained nevertheless because it 
is basic to all further development, and because Helfenstein’s paper is not easily accessible. 
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Fig. 4. Forcing amplitude factor B curves of simple elliptic response for a hard spring system. 
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Fic. 5. Forcing amplitude factor B curves of simple elliptic response for a soft spring system. 
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5 and those in Fig. 4 with k* < 0.5 cross the p-axis once. Near the points of crossing, 
such as the points J and J/, the value of B is extremely small. However, the response 
amplitude A is given by (14) and is finite. This means that near the crossing points 
extremely small external forcing amplitudes are capable of maintaining the response. 
This is the phenomenon of resonance. At any of the crossing points the external forcing 
required is zero. Here we have the case of free vibration. Thus, in non-linear systems 
such as those under discussion, the free vibrations emerge quite naturally as special 
cases of forced vibrations. 

3.2 Two-term elliptic forcing. Results similar to those obtained in Sec. 3.1 can 
be readily found for other cases. For example, let us require that the forcing function 


F(r) be B\x + B,x*. The equation of motion may then be put as, 
dx 3 / 
qt (1 — Bx + (6 — B,)x = 0, (20) 
aT 


and its solution can be easily found. It is 


x= Acn(pr,k), (21) 
where 
p = [1 — B, + (8 — B,)A’]'”, (21a) 
k = {(8 — B,)A?/2{1 — B, + (8 — B,)A’]}"”. (21b) 
The favored forcing function is 
F(r) = B,A en (pr, k) + B,A® en’ (pr, k). (22) 


Here we note that five parameters are involved in the solution, namely: A, B, , B; , p, 
and k, and there are two relationships between them. Therefore any three of the five, 
except the group of B, , p, and k, ean be arbitrarily chosen. A simple elliptic response like 
(21) is assured when the other two parameters are determined from (21a) and (21b). 
For the convenience of future applications we shall write out the favored elliptic forcing 
function, assuming that p, k and A have been chosen arbitrarily. This results in 


F(r A{l — p°(1 — 2k°)] en (pr, k) + A(BA® — 2prk’) en" (pr, k). (23) 


Another way of looking at the solution is the following: Let us denote B, A by C, and 
B,A®* by C,; . They are respectively the coefficients of the linear and cubic terms in (23). 
Combining C, and C; by eliminating A we obtain 


1/277 tte? BC? a 
°°, = —> Se ze SO — / 
7 “Ts - { — pl — 2k] 2pk } (24) 





as the necessary relationship between C, , C; , p, and k of the forcing function in order 
to have a simple elliptic response. 

It should be pointed out that although the exact solutions discussed here require a 
special class of forcing functions, that class itself is fairly general because of the freedom 
to choose any three of the parameters C, , C; , p, and k at will. 

3.3 Multiple-term elliptic response. The inverse method employed to find exact 
solutions in the last two sections can be extended to those consisting of more than one 
elliptic term. For example, if we require the solution to be 


a = A, en (pr, k) + A; en’ (pr, hb), (25) 
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the necessary forcing function can be found by substituting (25) into the left side of 
(11). Before doing this it is, however, convenient to perform a transformation on (11). 
Putting y = en (pr, k) and X(y) = 2z(r), Eq. (11) becomes 

eX 


pil — kh’) — (1 — 2k’)y’? — k’y*] 
dy (26) 


ees Dis weer 8 : = 
+ p[-(1 — 2h\y — 2k'y"] Ge + X + BX* = Gy), 


where G(y) = F(r). 
In our case (25), X = A,y + A,y*. When this is substituted in (26) and all terms 
are grouped in powers of y we get 


= Dy + Dy? + Dsy’? + Dry’ + Day’, (27) 


Gly 
where 

D, = [1 — p’(1 — 2k*))A, + 6p'°(1 — Az, 

D; = [@Ai — 2p’k*]A, + [1 — 9p°(1 — 2k°)]A; , 

D; = 3[BA; — 4p k’JA; , 

D; = 3BA,A; ; 


D, = BA;. 

Again, in this solution four parameters out of a total of nine may be arbitrarily 
chosen. An application of this exact solution will be discussed in the next section. It is 
evident how this result can easily be extended to cases where X(y) is of power higher 
than three. In general, if X(y) is of power 2n — 1 there will be 4n + 1 parameters and 
3n — 1 relations between them. Hence n + 2 parameters may be chosen arbitrarily. 

4, Simple harmonic forcing. Up to now we have only discussed some exact solutions 
under external forcing consisting of favored elliptic functions. When the forcing function 


is simple harmonic the equation of motion is 


dz . 3 2r 
— +2+ Br = P, cos T r = Py, coswr, (28) 
aT 


where 7’ is the forcing period, and w the forcing frequency. In this instance an exact 
solution in a simple form is hardly to be expected. Various procedures are available to 
find the approximate steady state solution. 

One method consists in assuming as an approximate solution a suitable function 
with some undetermined parameters, and then defining an error function e(r) as, 


2(r) 


dz ; 2r 
e(r) = oa +2-+ B2 — P, cos mT (29) 
dr i 


and, finally, determining the parameters such that the total error specified in some 
fashion be a minimum. From what we have learned in the preceding sections it seems 
natural to take a suitable elliptic function as the approximate solution. Let us put 
z(r) = A cn (pr, k). The error function is then given by 


; 2r 
e(r) = AB en (pr, k) — Po cos 77 (30) 
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provided that Eqs. (14) are satisfied. Furthermore, to assure the same period for the 
response and the forcing function we must have 


p = 4K(k)/T. (31) 


Here we have four parameters A, B, p, and k, but there are three conditions, Eqs. 

14) and Eq. (81), to be satisfied. Only one of them is completely at our disposal. This 

one can be determined by using any of the several minimum total error criteria, such as 

the colocation method, the Galerkin’s method, or the least square method. Solutions 

obtainable in this fashion were computed. When k’ is small, which means small | 3A? |, 
these methods yield the well-known result 


Po 


w = 1+ 38A? — r 


When a multiple-term elliptic response is used for the approximate solution 2(r), 
proper formulation of the minimum total error can be made in the usual manner without 
difficulty, although the actual calculation becomes very involved, as usual. 

We shall not reproduce these solutions here. Instead we wish to derive the approximate 
solutions of various orders of the equation (28) entirely from the various exact solutions 
obtained in Sec. 3. The central idea is to compare the orders of magnitude of the terms 
appearing in the exact solutions. It will be shown that approximate solutions of different 
orders will result quite naturally from this kind of analysis. 

4.1 First order approximation. Irom Sec. 3.1 it is seen that the response will be 

A en (pr, k), provided that the forcing function is of the form of (16). This favored 
elliptic forcing function will be rewritten here as F,(r). 


Fy(r) = A[l — p*(1 — 2k°)] en (pr, bk). (32) 
The function en may be expanded into the Fourier series 


‘ ‘ n+1/2 
(2n + l)apr 2r @q 


en (pr. / Re Qon+1 COS — us ji. * SS See 
De ont 2K , ome” 6S we 





(33) 


In this expansion it may be easily verified that the magnitudes of the various components 


are of the orders: 1, k’, k*, --- when k° is small in comparison with unity. Hence, outside 
of a common factor A the orders of magnitude of the leading parts of various harmonic 
components in /,(r) are respectively: (1 — p*), (1 — p*)k’, (1 — p*)k*, «++ . From 


14), 8BA~ is known to be of the order k’. 
Next, let us see what error we commit if we replace the simple harmonic forcing 
function P, cos wr by F,(r). When 1 — p’ is restricted to be of the order of k”, the first 


harmonic in /’,(r) will be of the order k° and the rest will all be of higher order. By setting 
the coefficient of the first harmonic equal to Py , we find 
Afl —p(l - 2k”) ] = P,. (34) 


We insure an accuracy of k° if we replace P, cos wr by F,(r), provided of course that k* 
is small compared with unity. In order to have the same periodicity for the response and 
the forcing we must have 

p = WwK/r. (35) 


Equations (34), (35) and the first equation of (14) completely determine the values of 
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A, p, and k, thus the first order approximate solution, for any set of assigned values of 
8, w, P, . Eliminating p from these three equations and expanding K in powers of k’, 


we get 


BA® — 2k*w*[1 + 3k? + O(k')] = 0, (36) 
' so Po * 
(l—w) — =FK'(1 — 250°) + O(k') = —- (37) 
- A 
Since, in this analysis, we limit 1 — p° to the order of k*, we find that 1 — w’ is also of 


the order k* and may be taken as ak” in a comparison of orders of magnitude, a being a 
quantity of the order unity. As the forcing function is approximated to the accuracy of 
k*, only terms up to that order will be retained in (36) and (37). These equations now 


become 


BA* — 2k* = 0, (36a) 
Q 2 
+ ee oe a7 
ak i 5) a A ° (37a) 
From them we find 
:.. Pe 2 Po 
»o =1+ $8 A° ~ A », or, w=1t+ 38 A° ~ A . (38) 
The response is x = A en (pr, k). Again, retaining only terms up to k° in its Fourier 


expansion, one gets 


l 
t= Po + 39 


BA*(—coswr + cos aun). (39) 

Equations (38) and (39) are exactly those obtained from the perturbation method, 
[8]. Although this first order solution result is not new, the method through which it is 
obtained here is novel in that it is obtained as a special case from an exact solution of 
forced oscillation. Throughout the derivation one thing stands out. In the perturbation 
method the perturbation parameter usually used is the coefficient 8 of the non-linear 
term. The present analysis indicates that the controlling parameter is k* or 8A*/2 rather 
than 8. Physically, this is to be expected. Even when 2 is not small, but A is restricted to 
to be small, the first order solution is still valid. 

Because of the requirement that 1 — w* be small, the solution obtained is only valid 
in the neighborhood of linear free vibrations. Let us now drop this requirement and 
assume that 1 — p’ may be of the order unity. In this case it is necessary to start with 
the exact solution found in Sec. 3.2. The response is again x = A en (pr, k) and the 
favored elliptic forcing function is given by (23) and will be denoted here as F,(r) 

The function cn* may be expanded into a Fourier series. 


(¢ 
-> be... cog ott _lepr 


en’(pr, k ; , 
— ans (40) 
7” (z)\ 
- = “RE Dk? — 2 meee 
Seas oe 2 1 + (2n + 1)° oR ber rs ot 


When k’ is small the magnitudes of the various components are: 1, 1, k’, k*, --- . It will 
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be shown later a posteriori that when k’ is small, 8A” is of the order k’. Keeping this in 
mind we find that, apart from a common factor A, the orders of magnitude of various 
harmonic components in F,(r) are: 

(1 — p’), i—p)k, (1—p)k, -++ from en term. 
yr. rr a k* : -++ from en* term, 


If we again neglect terms of order k* or higher we have two harmonics from the en and 
two from the en’. By reasoning similar to that used previously we can approximate the 
foreing function Py cos wr by F,(r) with an accuracy up to k’, if we set 


[1 — p(l — 2k’)Ja, + (@A* — 2k’*p*)b, = P./A, (41) 
[1 — p(1 — 2k*)Ja; + (@A* — 2k’p’)b, = 0. 


(nother necessary condition is that concerning the periodicity given by (35). Equations 
(41) and (35) completely determine the approximate solution, provided that k thus 
found satisfies the order-of-magnitude requirement. 

Substituting the coefficients a, , a; . b, , and 6, from (33) and (40) into (41), eliminating 
expanding K into powers of k*, and retaining only terms up to k’, we get 


> 


p by using (35), 
l1—w) — (1 — 25w°)k'/16 + 3(BA° — Qk’w)/4 = P,/A, (42) 


l wk’ /4 + (BA° — Qo) = 0. 


Eliminating i” from these two, we have 


’ Po , (jw — 1 
l a eS a m BA ieee“ = 0. (43) 

A %" — | 
After A has been determined from (43) for given 8, w, and P, , k” may be found from (42). 
k° = 48A*/(Qe" — 1). (44) 

The response is given by the following expression. 
x a me 
x = Ascoswr + — (—coswr + Cos 3w7r) ?- (45) 
16 J 

Near w = 1/3 the solution is probably not valid because of the large value of k° which 
invalidates the analysis. This solution also requires that | 6A* | be small in comparison 
with unity. It follows from (43) that 1 — w* — P,/A has to be small, implying small 


deviations from the linear forced response. The solution obviously represents a per- 
turbation from the linear forced vibration. It is interesting to note that in most of the 
frequency ranges (43) does not differ very much from (38). 

4.2 Higher order approximations. In the exact solution under two-term favored 
elliptic forcing we have three parameters completely at our disposal. However, in approxi- 
mating a simple harmonic forcing by such an elliptic forcing function, one relation 
between these three is necessary on account of the periodicity requirement. With only 
two parameters at our disposal we can not hope to obtain any approximation beyond 
the zeroth order and the order of k°. In order to get approximations of higher order we 
resort to the exact solutions of multiple-term elliptic response discussed in Sec. 3.3. For 
example, for the second order approximation we may use the response given by (25). 
The necessary forcing function is rewritten here as F;(r). 
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F3(r) = D, en + D, en*® + D; en® + D; en’ + D, cn’, (46) 


where en stands for cn (pr, k) and the D’s are given by (27). 

Let us assume that 6A’ is of the order k’ and A; is small in comparison with A, so 
that (A;/A,) is also of the order k’. We shall denote (A;/A;) by a;k’ in order-of-magnitude 
calculations, a, being a quantity of order unity. Under these circumstances the co- 
efficients D, and D, will be of the orders k°A and k*A respectively. If we retain only 
terms with powers of up to k*, the D; and D, terms need not be considered. The orders 
of D, , D; , and D; are respectively A, k°A and k*A. 

The function cen® may be expanded into a Fourier series. Let us denote the Fourier 
coefficients by ¢2,4; . They are respectively of the order 1, 1, 1, k’, --- when k? is small. 
Again by comparing the orders of magnitude of various terms in F;(r) we can approxi- 
mate the given forcing function P, cos wr by F;(r) to the accuracy of k*, provided that 
we have 

D,a, + D3b, + Dc, = P, ) 
D,a; + D;b; + Dc; = 0, (47) 
D,a; + D;b; + Des = 0. 


Retaining only the necessary powers of k in a’s, b’s and c’s to maintain a uniform accuracy 
of k* in (47), we have 


= pote 2p 4 OM 
a= =” i” 
a= +k + O18) 
‘ta Tae" Ok), 
_ . 
a,=sk + Ok), 
200 


bh =3-— Sy 4 O(k'), 


32 
1 3 - | x 4 
b, = 7 T 64 k + O(k), (48) 
) 
b, = 2 k? + O(k*'), 
64 
( ~ 
¢ = 0 + O(k), 
16 
ce = ~ + 0(k), 
16 : 


I] ae 
Cc; = 16 + O(k*). 


Substituting these coefficients and the expressions for D’s from (27) into (47) and retain- 
ing only terms up to k*, we get 
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1 —w* + 3847 + - (1 — w*)(12a3 — 1) + Pai ‘BAi(2003 — 1) 


_ soe k (1 => @ *) (Sas + 3) = P./A, ’ 
(49) 


BA; + 4k°(1 — 9w")(1 + 4a) +2 - ak*BA? + 7 k*(1 — 9w°)(2 + 3a;) = 0, 


K'BAU(L + das) + 75 kY(L — 25u")(1 + 12a) = 0. 


Here the periodicity requirement (35) has been used to eliminate p* in the D’s. These 
three equations completely determine the second order approximation. In principle at 
least, the values of A, , a, , and k” can be determined for any set of assigned values of 
P, , w, and 8. The solution is valid only when the resulting k’ is small and a; is of the 
order unity. Once the values of A, , a; and k’ are known, the response is given with an 


accuracy up to k* by 


f fy 4 
£ aqf) a 2( — 16 — ta,] a k (- 256 32% a.) | COS WT +. E (4 + tas) 
+ k (4 4 3 )| > s 3. ’ + | e( ] 3 )| 20S = ‘A. Ok” . 
39 = 64 a3 COS bWT * 556 64 Qs COS OWT (K) 


lor approximations of higher order than the second, similar analyses may be made 
by employing solutions with three or more terms in the response. No essential difficulties 
will arise except the complexity of the mathematical computation. 

Before concluding this section we summarize the essential features of the analysis 
presented here. Our aim is to find a steady state solution of (28). The general approach 
given here is, first to find some exact solutions of slightly different problems; in the 
present problem, namely: elliptic responses under favored elliptic forcings. At the first 
sight one may think that these exact solutions are very restricted. In actuality they 
represent a very wide class of solutions because of the freedom to assign arbitrary values 
to the parameters available at our disposal. Indeed, approximate solutions of various 
orders of (28) are but special cases of these exact solutions when k’ is small. Furthermore 
they can be obtained simply by comparing the orders of magnitude of various terms in 
the expansions without having to solve new differential equations which appear in the 
perturbation method. 

It is believed that this type of approach may find application in solving other non- 

linear problems. 
5. Subharmonic solutions and simple elliptic subharmonics. Consider again the simple 
elliptic response solution under two-term elliptic forcing, discussed in Sec. 3.2. The 
functions cn and cn* may be expanded into Fourier series according to (33) and (40). 
The lowest components are respectively: 


(50) 





2r q'” Tp 
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If we set these two equal in magnitude and opposite in sign, the fundamental component 


vanishes and there results 


2a, oes. 2. (5 ) : 
BA okt 21? \OK (51) 
Equation (51) expressed in terms of p, k, and A implies 
as = e , > > T i ~ 
BA* = 2k’p’ — 2k’*[1 — pil — 2k)] | 2%" —1+ (=t-) | (52) 


Under this condition the forcing function is devoid of the fundamental component and 
yet the same component in the response is predominant. This is the well-known sub- 


harmonic phenomenon. 
When k = 0, we have 


B, 1l—p, B, = —BA* ==s(1—-p), ait 
o { ed) 


x = A Cos pr, F(r) = —}3A(1 — p’) cos 3pr, 
which is the case of simple subharmonic of order 1/3. 
Similar analyses can be made for solutions with the sn function. It is to be noted 


here that although (52) leads to subharmonic solutions which are neither “pure” nor 


“simple”’ as defined by Rosenberg [6], they are nevertheless exact solutions. 

5.1 Simple elliptic subharmonic solutions. Analogous to the simple subharmonics 
we can define simple elliptic subharmonic solution of order 1/r as that x = A cn (pr, k) 
excited by a forcing function F(r) = C en (rpr, k). The condition under which such 
solutions exist can again be obtained by either the inverse method [7] or the technique 
of transformation described in [4]. 

Let us call cn (v, k) as z and define a function H(z) according to 


H,(z) = en (nv, k). (54) 
This is a generalization of the Tchebycheff’s polynomials of the first kind. The general 
function H(z) may be found by the addition theorems of elliptic functions. 


Next let us consider a solution 7 = cn (u/r, k). It satisfies the following differential 


equation 
d° ] Qk ' ar 
1 + >I — 2k°")n + = & 0 (55) 
duo ar Tr” 
If we add to this a trivial identity 
gH ,(n) g en (u, k) (56) 
we get a differential equation 
d ] . 2k ao 
— + -5(1 — 2k°)n + —> 1 + GH,(n) g en (u, k) (57) 
OU a r 


which represents a special non-linear forced oscillation under an elliptic forcing. The 
response cn (u/r, k) is also a simple elliptic function but of a period r times as large as 
the forcing function. This is defined as a simple elliptic subharmonic of order 1/r. 


In general H(z) will not be a finite polynomial but rather a more general rational 
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function in z. For example, H,(z) is given by 


2[1 — 48° + k*(6s") + k*(—48° + 8°)] 
1 + k*(4s° — 6s*) + k*(48° — 38°)’ 





H,(z) = (58) 
where s” sn” (v, k) = 1 — 2. When the modulus k° is small, H,,(z) can however be 
expanded into a series in powers of k*. Thus, for n = 3 we have 


y oo we Pr 2,2 on so" = 92" 42° 
H(z 2[(—3 + 42 ik°2"(—1 + 6 + 4z) (59) 
+ 4k*?(1 — 27)°(—1 + 327 + 122° — 162°) + O(k’)]. 


Substituting this into (57) with r = 3 and retaining only terms up to k’, we get 
eS © 


d°n 1-2, ) mt a a ee 


du \) 
— 36k’gn’ + 16k’gn° = g en (u, h), 


as the differential equation which possesses a simple elliptic subharmonic solution of 
the order 1/3. The solution is of the accuracy k’, where k’ is assumed to be small. 

When & 0, the function H,,(z) reduces to the usual Tchebycheff’s polynomials of 
the first kind and Eq. (57) becomes the one given by Rosenberg in [6]. 
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ny / 


(Continued from p. 37 


solved by the intelligent application of simple arithmetic’’. This policy has led to a choice of problems 
which require only a direct application of the discussed methods and formulas, and which therefore 
serve as useful exercises in probability algebra. Unfortunately, there is a corresponding artificiality to 
the problems, and there are very few which are either of a really imaginative nature or which lead to 
thought concerning their deduction from a real situation. On the other hand, this feature would not 
detract too much from the use of the book as a companion text in a ‘‘case-work”’ kind of course. 
The first chapter provides a survey o f probability concepts; although the authors do not hesitate 
] 
i 


to discuss “sample spaces’ or “interaction and union of derived events’’, they rather paradoxically omit 


the deriva 
dark comments as to the inherent mathematical difficulty of probability do not seem pedagogically 
useful. the chapter on Allocation, the simplex method is presented rather clearly, but with some 
mathematical awkwardness (cf. the incorrect statement of the key theorem of linear programming on 
p. 226.) As a final minor objection, the point of the “principle of optimality” (p. 272) in the chapter on 


tion of the Poisson distribution formula to be used dureashout the book. Also, the occasional 


dynamic programming is quite lost on the reviewer. 
Cart E,. PEARSON 


Instationdre Wdrmespannungen. By Heinz Parkus. Springer-Verlag, Vienna, 1959. 


- H. €0 0% 
165 pp. 39 Od. 


This book is a sequel to a previous volume on steady-state thermal stresses by E. Melan and the 
tise by providing a fairly extensive and up-to-date account 
of available inves tigati ms of thermal stresses and deformations due to time-dependent temperature 
fields. While most of the material selected here vere to classical elasticity theory, the book ends with 
an excursion into the fields of viscoelastic and elasto-plastic solids. 

The opening chapter contains certain elements of the linear theory of thermoelasticity for homogene- 


present author. It supplements the earlier trea 


ous, isotropic media. As in the preceding volume, however, the primary emphasis is not so much on a sys- 
tematic and complete exposition of basic theory as on a careful and detailed discussion of solutions to 
particular problems. Chapter II deals with investigations of transient temperature stresses pertaining 
primarily to an unbounded domain, the half-space, the sphere, and the circular cylinder. Chapter III 
effect of periodic temperature changes, whereas problems involving moving heat 
Chapter IV. These quasi-static analyses are followed by a chapter on inertia 
Thermal 


is devoted to the 
sources are treated ir 
effects in transient thermoelasticity which presents a resumé of recent studies in this area. 
stresses in viscoelastic and elasto-plastic media are taken up in two concluding chapters. 

Some of the results included by the author are in the form of improper integrals. It is difficult to 
appreciate the usefulness and significance of such solutions to boundary-value problems in the absence 
of adequate numerical evaluations which are missing in several important instances. On the other hand, 
integral representations of singular solutions which correspond to thermoelastic nuclei are bound to be 
of limited interest in the further development of thermoelasticity theory. 

A brief discussion of thermoelastic coupling effects might well have been appropriate in an enter- 
prise of this kind. Further, the chapters on temperature stresses in elastic materials are rather sketchy 
and reflect the ritualistic character of our present approach to such problems—an approach that is 
for the most part safely removed from the complexities of the actual physical situation at which it aims. 
A more realistic treatment of thermo-inelasticity cannot escape the essential nonlinearities attached to 
this topic and must come to grips with the temperature-dependence of the physical parameters which 
govern the thermal and mechanical behavior of such materials at elevated temperatures. 

It is hardly fair, however, to blame the author, who has made a very welcome contribution to a very 
timely subject, for inadequacies in the present state of this subject. The book is clearly written, care- 
fully edited, and admirably printed. The comprehensive bibliography appearing at the end of the volume 
will be greatly appreciated by every interested reader 


E11 STERNBERG 
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TORSION AND EXTENSION OF HELICOIDAL SHELLS* 


JAMES K. KNOWLES (California Institute of Technology) 
AND 
ERIC REISSNER (Massachusetts Institute of Technology) 


1. Introduction. The present paper is concerned with the problem of rotationally 
symmetric deformations of thin elastic shells the middle surface of which is a portion of 
a right helicoid. Particular consideration is given to the problem of a uniformly pre- 
twisted thin strip which is acted upon by tractions which result in equal and opposite 
axial forces F, and equal and opposite axial torques 7’ (Fig. 1). 

The differential equations for stresses and deformations of thin homogeneous, isotropic 
helicoidal shells, as used here, have been derived elsewhere [1]. In what follows they are 
employed in the form which they assume for rotationally symmetric states of stress and 
strain. Insofar as the problem of the pretwisted strip is concerned one of the essential 
aspects of the analysis is the connection between rotationally symmetric states of strain 
which depend on states of displacement which are not rotationally symmetric. The details 
of this connection are established in the present paper. 

Of particular interest in the problem of the pretwisted strip are the relations between 
the applied force and torque on the one hand and the angle of elastic twist and the 
relative axial extension on the other hand. In this connection certain explicit results 
are presented which generalize earlier work of Chen Chu [2] regarding the torsional 
rigidity of the strip. 

2. Equations for helicoidal shells. Let r, 6, z be cylindrical coordinates and let 


z= age (2.1) 


be the equation of the middle surface of the shell. The constant 27a is the pitch of the 
helicoidal middle surface. The parametric curves r = constant and @ = constant on the 
middle surface of the shell form an orthogonal net but are not the lines of curvature. 
The state of stress in the helicoidal shell (2.1) is described by stress resultants N, , 
No, Neo, No , Q, and Q, and stress couples M, , M, , M,, and M,, referred to tangential 
and normal directions at the edges of an element of the shell (Tig. 2). The differential 


equations of equilibrium of an element of the shell are [1], 





0 ON " ra a ‘ 
(aN,) + —* —-N,+-Q, = 0, (2.2) 
} 00 Y 
re) ; ON, = a 
— (aN,») + — +-N, +-Q, = 0, (2.3) 
or 00 a a 
ra] , oc Jo a a ov ¢ 
. (aQ),) 4. = = (N ro +: NA 6,) = 0, (2.4) 
or 00 a 
re] OM,, s ” 
—(aM,) + — — —-M, — aQ, = 0, (2.5) 
or 06 a 
*Received August 4, 1958; typographically revised manuscript received December 29, 1958. A report 
on work carried out under the terms of an Office of Naval Research Contract at the Massachusetts 


Institute of Technology. 
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Fic. 1 

















Fic. 2 Element of helicoidal! shell 


9 : IM , r 
: (aM 4) + om a i Mo, — aQ, =,0, (2.6) 
Or og Qa 
Ny» — Ne +3 (Ms — M,) =,0. (2.7) 
a 


The quantity a is defined by 
a = (a +r)”. (2.8) 





—— 
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The system (2.2) to (2.7) is completed by a system of stress strain relations which 
here is taken in a form corresponding to the relations of Fliigge [3] and Byrne [4] for 
lines of curvature coordinates and which is [1], 


N, = 3 ee, (, + ve) — a5 Eh | + -(-)5@- «|, (2.9) 
Ny =; a (@ +) — 4 ae = | -(l-) 4 - «| , (2.10) 
Nve a <0 — me (1 + xt + (3 — yf], (2.11) 

No SAY + ¥e0 — = sea Ia (1 + va + (3 — ») xP], (2.12) 

M, == re > (<r + nt — 25254, ) (2.13) 

Mo = TH on (xs + ot -* —rs he ), (2.14) 

Mw = zy En . I 5 Tt — a (eo + 6, |, (2.15) 

Me = 7 cu = E = dy + ve |: (2.16) 


The strain quantities ¢, y, x and 7 are expressed in terms of radial, circumferential 
and axial components of middle surface displacement u, v and w as follows; 


Ou oy 
6.=—, (2.17) 
or 
r ov a ow r 
= 3 — + a + 8 (2.18) 
a 06 a 06 a , 
1 Ou r ov aow v 
vent gt Sy tS. (2.19) 
a 00 a Or a or a 
r ow a’ ow adv ar dv a a ou 
x* ont oe a -_ + re ees Oe + Te =e (2.20) 
a or a or a or a’ or a a’ 06 
is Fil r dw ada v = 1 dw a du (2.21) 
a oo a oe a or a’ 00 
" 2r aw 2a dv 2(r? — a’) dw far dv dar (2.29) 
- -= — — = —— — -;- - = tl. 2.46 
a oro@ a oroé a’ 06 a’ ae 


lor the formulation of boundary conditions it is necessary to have expressions for 
effective edge stress resultants which include the action of the edge twisting moments 
and their derivatives, insofar as they are statically equivalent to distributions of edge 
forces. Radial, circumferential and axial components of these effective edge stress re- 
sultants are, along edges r = constant, 


R,=N,+- Mw, (2.23) 
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Fo a a oM,, 
H, =-N,, —-Q, -3a— (2.24) 
a a a 006 
, r r OM,, . 
Z,=-N.+-Q, +3; 2.25) 
a a a 00 
and, along edges 6 = constant, 
; a 
R, = Noe t+—aM,, , 
T (2.26) 
/ \ 
7 a Oo a = 
H, = -N,-—-Q, —- Uo.) , 2.27) 
a a OT \a@ / 
a a 7 Qo (rT 
Zo =-Not-Q+ ( ar,,) 2.28 
a a OT \@ 
In addition this procedure introduces concentrated corner forces of magnitude 
+ (M,, + M,,) in the direction of the normal to the middle surface of the shell, with 
components + (7r/a) (M,, + M/;,) and + (a/a) (M,, + M,,) in the axial and circumfer- 


ential directions, respectively. 
3. Boundary conditions for the problem of the pretwisted strip. We consider a 


helicoidal shell with edge coordinates r = + band 6 = + @,. The usual polar coordinate 
interpretation is attached to the meaning of negative values of r; the point (— r, 4) is 
the image of the point (r, 6) under reflection in the z-axis. 

We assume that the edges r = + b are free of stress and that the edges @ + 9 


are acted upon by forces F and torques 7’ in accordance with Fig. 1. We then have the 


following system of boundary conditions along edges r = constant: 
r= +b: k.=H M, = ®, 3.1) 


Along edges 6 = constant the boundary conditions are taken in a form which insures 
that a rotationally symmetric state of stress will exist in the strip. Thus we prescribe, 
atdé = + 6 


[ de a f . oe ue) | F, 3.2) 
[Hy dr + k M.,+M | a’ 3 3) 
| N,. dr =0 5 4) 


3.6) 


—, 
— 
| 
S ~ 
= 
+ in 
Ss 
— 


4. Displacements and stresses for the pretwisted strip. Considerations of symme- 
try indicate that displacements for the pretwisted strip problem should be of the form 
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u = u(r), v = v(r)6 and w = wp(r)é. The requirement that the components of strain 
(2.17) to (2.22) be independent of @ leads to the conclusion that admissible displace- 
ments are of the form 
u = u(r), v = kré, w= k,6, (4.1) 

where k, and k, are constants. 

lor displacements of the form (4.1) we have further that «* = x} = y,, = 0. In 
view of Eqs. (2.11) to (2.14) this means that the displacements (4.1) are associated with 
the vanishing of four resultants and couples, 


N= N,, =M, = M,=0. (4.2) 


The equilibrium equation (2.3) then implies 
Q, = 0, (4.3) 


c 


g system of differential equations for VN, , Ng , M., , Mo, , Qe , and u 


., = , 
and the followin 


remains 


alaN,)’ — rN, + aQ, = 0, (4.4) 
alaM,,)’ + rM,, — a°'Q, = 0, t.5) 
] a Ek} 1 
\ € ve sir — (1 — = te — €) t.6) 
= a 12(1 — vr’) a 
ke) a Eh* a 
\ (€g + ve.) — a | 7* — (1 — v) -s (ey — e,) 1.7) 
| I a 1211 — v) a 
Eh l—vp a 
1 = -|- r* — — (e, €,) (4.8 
q 1 — >) 9 7 (eg + ve,) |, {.8) 
Eh’ l—p a 
Me, 2 mange | Se gt oe Sale, 4 pa) b (4.9 
] 120 — 2) 9 T AC + ve .9) 
In these equations primes denote differentiation with respect to r, and 
r r’ a 
‘. =e, €, sutkhitkh-s, (4.10) 
a a Qa 
lar 2a(a® — r’) 2(r? — a’) 
m= uth t+ kh, >. (4.11) 
a Qa Qa 


Three of the four boundary conditions (3.1) are automatically satisfied and the fourth 
may be written in the form 
r=+b: aN,+aM,, = 0. (4.12) 


Of the six boundary conditions (3.2) to (3.7) four are automatically satisfied and the 
remaining two may be simplified (upon elimination of Q,) by suitable integration by 
parts to read 


a Qa 


r= | (: Net 5" Ma +5 M,,) dr, (4.13) 


r=f (x +°S% M445 Mu) ar (4.14) 
© b Qa 


a Q 
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5. Non-dimensionalization and simplification. In Eq. (4.1) we may write 


ky = do, Ka = 06, (5.1) 


where w and 6 are respectively the angle of twist and the axial extension, both per unit 


of axial length. 
We further introduce a dimensionless displacement and dimensionless resultants 


and couples as follows: 


Uo = ; (5.2) 
_ N, _ Nz a (5 2) 
Me Eh? Bh? = BR? _ 
My _ Ms, > 

me = Eh? MN, = Eh? 5. 


The quantities up , n, g, and m are considered as functions of a dimensionless coordinate 


p defined by 


i a ~ 
) = = (0.0 
mb 
We set finally 
b — 
A =———e (2.0) 
a 


The parameter \ measures the pretwist of the strip and vanishes for an untwisted plate 
located in the xz-plane. 
Introduction of (5.1) to (5.6) into Eqs. (4.4) to (4.11) transforms these equations to 


the following form 


fa +2? dp ! 4 0 (5.7) 
“a | 7. 8 } 7 = } Ile . > > 572 96 = \o 
dp r ie (Il +X p)* saith (1 + dp)’ ~ l : 
a. ahi x Haye ss 
ae AE ee ig “M,¢] - —— 2 F172 Me, — (1 + Np) a = (), (9.3) 
ad (1 + Xp) 


‘ l h* r (l—vp)rA, 
— yp )n : { Jeg “¢ 7 an br* oom nivale os = a) a (5.9 
| v)n, €, + ve ite! I + vp? € | (5.9) 


1 h d (1 — v)A . 
l—v)no = & +ve, — = = 7 ic - 3 (€ — € | , (5.10) 


12b° 1+ Xp" 1+ d°p 
12(1 — v)m,, = . > ~ by* = , er, (€g + ve,), (5.11) 
12(1 — v*)m,, = " = be* — ow (e, + ves), (5.12) 
where 
€é, = 7 ; = 7 rains Uo(p) + ean + rx, ; (5.13) 
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—4) p 1 — d*p Np? — 1 i 
se : BR Book Ff CS Be 5.12 
IT G+» oe 5 Uy + 2bw a] 2 oy 5 + 26 al Wa (5.14) 


The system (5.7) to (5.14) is to be solved subject to the boundary conditions (4.12), 
which take the dimensionless form 


p= cl: n, + - 7, ; x mo = 0. (5.15) 

The expressions (4.13) and (4.14) for the force F and the torque 7’ assume the form 

F = Ehb [ 7 r ee + ae ae : ae mia + 7 aa : a me dp, (5.16) 
T = Ehb’ [ | 3 Ma Ne + aie k, M,6 

+ a a" i h me | dp. (5.17) 


We now assume that the pretwist parameter ) is of order of magnitude unity and 
not large compared with unity, and that all dimensionless resultants and couples are of 
the same order of magnitude. Considering the fact that h’/b” < 1, we may then neglect 
certain of the terms in the system (5.7) to (5.17) and thus reduce it to the form 


Xp 


d 7 ) 7 . 
dp (1 + Xp)'’’n,] — Gj Teme = 0, (5.18) 
a 1 + d*p)'’?m,o] + dp me, = (1 + dX’)? G5 (5.19) 
dp ae? owe kee sili ; 
2 (5.20) 
l—yp 
ny = tte (5.21) 
l—vp 
— 2) =e oo Me ia 
m.« = 120 — %) | 2 br i+ Xp 5 (Ee + ve >|, (5.22) 
l—»,, d | a 
. =a !- — br* — ——- al, .23 
- eer | asd $y +0 one 
p= cl: n, = 0, (5.24) 
’ Ne » 
= Ehl ——*" 773 dp, .25 
I Enb | Cl + 2p de (5.25) 
Ap 1— 2) bh’ 
il = a va ¥ » \3/ 
1 ent [| ya eT Ne Be Me 
1 h? ® 
pe +n)” 3 3 me | dp. (5.26) 


6. Reduction of the differential equations. The resultants n, , n, , and the couples 
Mm,» and me, are expressed in terms of uo(p) by (5.20) to (5.23). The moment equilibrium 
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equation (5.19) serves to express g, in terms of Up : 


—vr dup 


1 + dp dp’ 


1211 — v’)q. = 


(4 — 3y)X"p dus 
(1+ A*p)° dp 
(1 — v)A‘p — 2(3 — vd (7 — 5v)d° ' 
+ bw — RE 5 (6.1) 
(1+Xp) (1 + Xp) 

The remaining condition (5.18) for radial force equilibrium provides the following 
differential equation for up . 
dup dp dup 
— x = ie ere wee i) 


de +1400 de tate 
== bof _ v)Ap = (d+ y)? | + 6 1 + y)X p (6.2) 





(1 +p) 


The boundary condition (5.24) then becomes 

duo vr" p vip v 

p= +1: _ 5-5 Up + bw — “3 6 = azz = 0. (6.3) 
dp 1+. °°  1+)'o 1+ ’p 


Thus u(p) is a solution of the boundary value problem (6.2) and (6.3); it will depend 


on w and 6. 
In view of the linearity of the boundary value problem we may write 


Uy = bwu, + buy (6.4) 
and split (6.2), (6.3) into the following two boundary value problems for u, and wz . 
(1 — v)Ap (1 + v)Ap sta 
hh = a ne (6.5) 
1+Xp l+2Xp 
du vr" p —vyr —_ 
p + | + ss Ul ay (06.0 
dp 1 -+2Xp 1+A 
(] + v)X\"p — 
Lu > 2 . 6.4) 
ray 
du vp -V 
p +1: oo 5 5 Us = Pts 6.3) 
dp l+Xp l+Ap 


where Z is the differential operator on the left side of (6.2). 

The differential equation Zu = 0 is, in different notation, the equation derived by 
Sanders [5] for helicoidal shell problems in which the displacements are independent 
of 6. It is possible to reduce this equation to hypergeometric form in a number of different 
ways, but no use will be made of this possibility in what follows. 

7. Influence coefficients. When the boundary value problems (6.5) to (6.8) have 
been solved, we shall have all quantities expressed in terms of known functions of p, 


and the constants w and 6. Thus by (5.20) to (5.23) and (6.4), 


, li > vip 
(l—y)n, = bo + - P. 3, + 5 
dp 1+ p 1+ Xp 


dus vr* p v | sy 
+ § 24 — —5 U renee | (7.4) 
| 1+ dp ” T 1+ Xp ; 
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du, dp Ap° 
(1 —v’)n, = bol » = + - 5 Uy, + | 


5 4+ Ne 2 
dus ’p 1 = 
+ fot Me t hs], (7.2) 











—yr du, (3 — 2)r*p i-)-—-@- one] 

o(1 — ome = bel ——A as _ B= 2" aes 
12(] y’)m, EF + 0 dp a+™ 2)? u, + (1 + Np) 

|  —vA_ du, (3 — 2y) n° p (1 — v)r\*p? — -3— ary = 

; + dp tue (1 + Xp")? + oS 

—r du, (2 — v)d\*p (l—»v) - Me] 

1 — v*)m, ot SS =F; 
12(1 v )mMe, bol 5 +} > dp (1 + n2p")? uU, + +xp) 








l+A\{p dp (1+ Po si (1 + 2°’) 


When these relations are introduced into the force and torque conditions (5.25), 
(5.26) we have two relations between the force F and torque 7 on the one hand and the 
angle of twist w and the axial extension 6 on the other hand. These relations may be 


+ | = duy rane = y)" - ot+ a st v)d*p f 4]. (7.4) 


written in the form 
F = (2Ehby,p;)6 + (ZEhb dyr.)w, (7.5) 
T = (2Ehb dy73)6 + (RED Ny 7. oe 2Gh* by)w, (7.6) 


where the dimensionless influence coefficients y are given by 


a 1 af y ¢ du a pus | ] 
ee  im-v ah Leer’ @& °° O+¥ar +e 


a! u + it dp, (7.8) 











l 
= l—vh, Cl aye (1 + \*p)*” (1+ Ap) 
[ du 7 p° p 
‘s a4 — ———3-n773 | d (7. 
v1 awa dp * G4 np)? 2 + G+ vp) , 7.9) 
du, dp Ap’ | “ 
4 mee +- — a an ee te 5 at De d ; (7.10 
_ UE Las 1+ x "i "dp | (+N + pyro” 
and finally 
[ _ dp 1 + (47/3) + (SX"/15) (7.11) 
4 (1 +-+xXe) Te Ss \é. 


In the integrals (7.7) to (7.11) use has been made of the fact that wu, and uw, are odd in p, 
so that the integrands are even and f!, = 2f). 
It will be shown in Sec. 8 that when \ = 0 we have yr; = ¥ru = ¥78 = Y¥ruo = Y = 1, 
so that (7.5) and (7.6) become 
F = 2Ehb, T = 3Gh’bo, 
corresponding to well-known results in the theory of extension and torsion of a flat plate 


by end loads. 
A straightforward application of Green’s formula, making use of the fact that u, 
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and w, are solutions of (6.5) to (6.8), shows that yr. = yrs , as is in fact required by the 


reciprocal work theorem of elasticity. 
In addition to the flexibility relations (7.5), (7.6), we have the associated inverse 

















equations 
= K.7T + KurF 
eae ie } (7.12) 
6 = K,,7T' oe K;rF 
It is readily shown from (7.5) and (7.6) that 
. 3 
Kur - 2Gh°t Diy(A ) + (4k b?/ 15Gh? wi on] 
a oe ee —)k(A) ; { an 
Kir = Kor = o6750) + GEb/15GR Nf) | [’ (7.13) 
ae __ksr() + | (3E b’/ 5Gh* xs "g(A) | 
oF ~ QEhbly(A) + (4EB°/15Gh)v7f() ] 
where the abbreviations 
Ja vy, — 2 Yeorrs 
f(A) =e 4 VT 4 Yr3 ’ 
k(d) = Y2e = Ue, (7.14) 
YF3 YFs 
keQ) =—, gd) = 77 
Yrs YF 
have been used. 
In the absence of axial forces, the torque 7’ and twist w are related by 
T = Iw, 
where the torsional rigidity J is derived from (7.12) and (7.13) as 
1 1 Hb 
= sa = 2G b / . } ) (F14e 
I K.; th? >a) +4 #Y x10) | (7.15) 


In the absence of torque the axial force F and the relative extension 6 are related by 
F = Ké, 
where the axial stiffness K is given by 


> A) — b* L5G / 
ee 1 — 2Ehb yA - (4E iO” Gh" AFA) 


Kor ks, r(A) 4 + (3Eb?/5Gh’) yr? g(r) (7.16) 


In the following section the first three terms in the power series expansions for 
yrs(\), Yro(A), Yrw(A), ¥(A), f(A), (A), Ksr(A), and g(A) are obtained by perturbation 
methods. 

8. Perturbation solutions. All quantities of physical interest have been expressed 
in terms of the solutions of the two boundary value problems (6.5) to (6.8). Inspection 
of these problems indicates that u,(p, A) is odd in \ and u,(p, A) is even in A, while both 
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are odd functions of p. For sufficiently small , we obtain solutions in the form 
uy(p, ) = duy"’(p) + Nur?'(p) + - °° ,| 
Us(p, A) = Us (p) + dus”? (p) fp see, J 


Introduction of these assumpticns into the boundary value problems (6.5) to (6.8) 
leads to a sequence of boundary value problems for the u;’)(p). These may be solved 
successively by repeated integration. The results of such calculations are 


(8.1) 








v l—y l v ye 5 1,3 
mh seat | ; ** (B+%+%) h 
fs . 9+ Sy (1-1 +9) 5 _ (135 + 383» + Sty + v) zs 
L 35 C«S 24 p 2520 
L OO%. (8.2) 
l-—v (1 /)° r 3 
—vp + E — — pt+~ ~ a aS + E (1 +») —v)p 
2 6 8 
| , > U+r*O5 +7) 5 ha 
+ 51 — 1 + 9" — * +H05 +9 yh 
| (1 ++)°1 — G1- 3%) a -*)(1 +») 5 
L 144 , 16 
_ - vy) 1 + v)(15 + ») By (1 — v*)(515 + 60r + vr’) aN 
240 . 5040 p 
4: Of>°). (8.3) 
6 T T T T 
J 
/ 
SF 7 
/ 
/ 4 
First Approximation i Mh 
Chen Chu) . 
<i "Micah F 
4 
ae 
1, Ir a 






Second 
2 Approximation4 











| a a ee “ll 
Zeroth Approximation 
(St. Venant) 
Oo l | 1 
.@) 0.1 0.2 0.3 0.4 0.5 


d 


Fic. 3 Torsional rigidity ratio vs. pretwist parameter for y = 1/3 and 2b/h = 10 
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Introduction of the perturbation solutions (8.2) and (8.3) into the integrals (7.7) 
to (7.10) gives the power series expansions in A of yrs , Yr. = Yrs, and yr. ; (A) may be 
directly expanded from (7.11). These expansions are 


3+ 4)... 29 + 88y + 56” 











— a 4 6 
vr = 1 i s Mt + O(d°) 
9+ 8... , 161 + 304 152,” 
yee =n = 1 Py 4 SE at +00), 
| (8.4) 
7 45 + 20v ). , 531 + 504» + 232", : 
te <b ey + aa Mt + O(r’), 





1-1) + Sy + O(r’). 


II 


y 


These in turn provide the expansions of f(A), k(A), kse(A) and g(A) defined in (7.14). 














3—4y.. , 427 — 152 8 
fai Ses + On, 
92 ae 2 
k = 1 we 6 1 x? 4+ 42 + dy rv! + O(a’), 
, a (8.5) 
kp =1— <1 —ynv+ 5. Bis = — * + O(r'), 
12 — 4... 828 — 486y + 52r’ ., 6 
= ]-—-—-——)\ rN d°). 
9 21 + 2835 + 00).| 





In particular the torsional rigidity J of formula (7.15) takes the form 


_ _ a 63 ++) 
r= 1 (1 ry + a + 


Eb? ( 33 — 4v., , 427 — 1527+ 8, )| 
C6} 86a 840 pitied ft 


| 





v7 


on 


where 
Ty = Theo = 3Gh*O (8.7) 


is the St. Venant torsional rigidity of a flat plate of thickness h and width 2b. If only 
the first term in each power series in parenthesis in (8.6) is retained, we obtain the Chen 
Chu approximation [2] 


2 
r=1(1+42%) (8.8) 


indicating the increase in torsional stiffness for small pretwist. Figure 3 compares the 
Chu approximation (8.8) with the second approximation (retaining \” terms in the two 
power series) and the third approximation (retaining \* terms). 

Corresponding results for the axial stiffness K of (7.16) are obtained by inserting 
the power series for y, f, ks and g into (7.16). There follows 
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, - »¢ iy 2 —_ 
K _ ( 7324 684 .-+) a SEB v(t _ 33-4), 
6 42 








Re 40) 15G h 
427 — 152y + 8’ ,, --)]{( _ 2-27), , 45 - 38%—-+",, ++) 
+ 840 + 3; *r 45 * 
3Eb..(, 12—4y., , 828 — 486+ 52”, yr 
+ 5G Px 7; ** 2835 + , (8.9) 
where 
Ko = Kyo0 = 2Ehb (8.10) 


is the axial stiffness of a flat plate according to plane stress. If only the first term is 
retained in each of the four power series in parentheses in (8.9) there follows what may 
be considered as an analogue of the Chu approximation; 


1 + (4E0°/15Gh?)X? 
1 + BEW/5Gh)X 


Figure 4 compares the first approximation (8.11) with the second and third approxi- 
mations obtained by retaining the d* and \‘ terms, respectively, in the four power series 
in parentheses in (8.9). 

In addition to the influence coefficients, the stress resultants and couples themselves 
may be calculated by introducing the expansions (8.2) and (8.3) into the expression 





K=>K, (8.11) 
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Fic. 4 Extensional stiffness ratio vs. pretwist parameter for »y = 1/3 and 2b/h = 10 
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(7.1) to (7.4) for n, , Np , Mm,» , and ms, . We obtain 


a 9+ 5 L-s 5 +P? ahs 
n, = bol —40 — p') + (2+ 7 - Ha 2*)a + | 
i s| —10 = p°)d? 3 (A+ + i= ba p° = 2 t 5) 4. | . (8.12) 


ive x a oe Oy + 5r° 2-—yv-yYr , 

~= be a ae (s i Sn aa ) 3 ( - oe le a 2 
Ng |» r + 12 p |X + 26 8 p 
378 + 57v +" 6\5 , , vy , 2+» 2\2 

+ SEES hes. | 4 fr - (G+ 25% ep 


38v +3 2-—vy—v ., 28 l7v + ¥* 
+ (2 oe SP ys ae ce tht +---], (8.13) 


- pty | 844+ 338v+0° , 
12(1 + »)m,» = bol — (4+ yp + (pte +34 one o')a aa | 


y+ O+%+» . é w + 3y 
+ -@torg (Ey WEete ley (REM TS 





ri Z 2 8 
6—yv—S8r-yr , 204 + 217y + 38° + v»” ,)\.; 
—— oe pjJA + ---1, 8.14) 


; 1+yp 57 + 5y , 
1211 + v)me, = bel 1 —3pr + fs. si = ? = p )n' + °° | 
+ 





1—2—3r , 101 + 82» + 5y . de 
2 —_ = p as ier we 5) + sof (8.15) 


If (7.12), (7.13) and (8.5) are used to express w and 6 in terms of F and T in (8.12) 
to (8.15), the expression for stress resultants and couples are obtained in terms of F 
and 7’. The corresponding expressions are lengthy and will not be given here. 
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WAVE PROPAGATION IN A COAXIAL SYSTEM* 
BY 
V. M. PAPADOPOULOS 


Brown University 


Abstract. A solution is obtained for the problem of the propagation of electro- 
magnetic waves in a semi-infinite flanged coaxial line with an infinite center conductor, 
in terms of an infinite set of coefficients which are determined by an infinite set of linear 
equations. The solution is discussed, in detail, in limiting cases which illustrate properties 
both of a thin vertical antenna on a plane perfectly conducting earth, and of a thick 
antenna fed by a low impedance line. Numerical results are given in these cases. The 
possibility of a solution for any excitation frequency is also discussed. 

Introduction. The use of an antenna for the purpose of radiating electromagnetic 
energy must involve the generation of this energy and its passage to the antenna from 
the generator along a transmission line. There is, however, very little theoretical work 
published in which an antenna is considered with its means of excitation. It might 
appear to be an advantage to examine the radiating system in isolation, and then to 
regard it as an impedance lumped at the end of a transmission line: the magnitude of 
such an impedance depends, however, on the transmission line parameters. It is therefore 
of interest to find whether there is any range of parameters for which the line and the 
antenna are substantially independent. The work done on this subject by King and others 
[3] implies that as we reduce the spacing of the transmission line to zero we may ex- 
trapolate from a series of experimental measurements of physical quantities to the limiting 
case of zero spacing; this value may be identified with theoretical results obtained by 
assuming, in the isolated radiating system, a delta-function excitation (sometimes 
called a “slice’’ generator) at the driving point of the antenna. This statement needs 
qualification, since the question of the existence of the physical quantity in the limit of 
zero spacing was not considered by King. We find that this limit does not indeed exist. 
Accordingly, in this paper we are concerned with an idealization of the cylindrical 
antenna problem in which the complete transmission circuit may be examined math- 
ematically. As an approach to the problem of the antenna of finite length, we shall 
consider a semi-infinite flanged coaxial line with its center conductor extended to an 
infinite length. A signal set up at one end of the line is partly reflected at the open end 
and partly radiated into the half-space outside the line. Since we are interested in 
antennae, we may describe the radiating system as an infinite vertical antenna of circular 
cross-section, standing on a horizontal, perfectly conducting ground of infinite extent, 
with the exciting signal fed in at the bottom. 

We shall further simplify the system by considering the case in which the field is 
radially symmetric about the axis of the cylinder, with no magnetic field component 
parallel to this axis. The field components are related by Maxwell’s equations; in this 
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problem all the non-zero field components may be written down in terms of the transverse 
component of the magnetic field tangential to the surface of the cyclinder, and this 
component satisfies the scalar wave equation. The normal derivatives of the magnetic 
field vanish on all the perfectly conducing surfaces. The field components in the half- 
space are related to the normal derivative of the magnetic field in the terminal plane of 
the coaxial line: this relationship may be found by a method in which cosine transforms 
are used. Since the field within the line may be expressed in terms of an infinite set of 
discrete modes (Marcuvitz [4]), it is possible by matching orthogonal components across 
the gap at the open end of the line to set up an infinite set of linear equations involving 
an infinite set of unknown coefficients. Each coefficient is related simply to the amplitude 
of the corresponding mode in the line; in turn, the set of coefficients determines all the 
field components completely and uniquely. To simplify the analysis, the free-space 
propagation constant is taken to have a small negative imaginary part which is later 
taken to be zero. The resulting solution is shown to satisfy the conditions of the problem. 

The solution of the infinite set of equations is simplified when the spacing of the 
coaxial line is small compared to the wavelength of the exciting signal. This simplifying 
condition is valid either when the outer diameter of the line is a small fraction of the 
wavelength or when the spacing between the conductors in the line is a small fraction 
of the line diameter. With this simplification, we find that the field within the line is 
very nearly that in an open-circuited line, and this approximate field leads us to a good 
value both for the admittance of the radiating system and for the energy density at 
large distances from the antenna. The effect of the simplification is to give us results 
not only for the thick cylindrical antenna with a small line spacing but also for the thin 
antenna with a physically plausible method of excitation. 

It should be emphasized that the results are appropriate to an idealized antenna of 
infinite length. The energy distribution may not be related to that of a finite antenna. 
On the other hand, the admittance of the infinite antenna is the limiting value of the 


admittance of long antennae. 
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Results associated with a similar geometry without a flange have been given by 
Marcuvitz [4]. It is implied that these results were obtained by a method involving the 
use of the Wiener-Hopf technique. In the present problem a similar technique may well 
be practical, but the author has not yet made a comparison. 

I. The geometry of the system to be considered is shown in Fig. 1. Cylindrical 
co-ordinates (r, @, z) are used. The axis of an infinite, circular, conducting cylinder of 
radius b > 0 is taken to be the z-axis. In the region z < 0 this cylinder is surrounded by 
a semi-infinite coaxial conducting surface of radius a (a > b), and this surface, r = a, 
z < 0, is terminated in the plane z = 0 by a perfectly conducting plane in the region 
r > a, z = 0. A time dependence exp (twé) is assumed throughout, and rationalized 
m.k.s. units are used. 

We are to consider the electromagnetic field in the region bounded by the conducting 
surfaces corresponding to a combination of axially symmetric transverse magnetic modes 
in the line region a > r > b, z < 0. This type of field is independent of the co-ordinate 
6, and it has the magnetic field components B, and B, both zero. From Maxwell’s equa- 
tions, it follows that the electric field components are related to the magnetic field 
component B, by the equations, 


7.2 
* p,=+2 @Bp, 
w ror 
ik? _ OB, (1) 
—E, = -—, 
w 0z 
E, == 0, 


where k is taken with a small negative imaginary part (k = k, — tk; ,k, , k, > 0), and 
B, must satisfy the equation, 
Set Pe , t) 
(2.+12+4+0 -3)p =o. (2) 
Since the tangential component of the electric field vanishes at a perfectly conducting 
surface, B, must satisfy the boundary conditions 


2 WB) =0 on r=b, (3) 
or 
0 
—(rB,.) =0 on r=a, z<0, (4) 
or 
0 
a, (Be) =0 on z=0, r>a. (5) 


We shall assume the total field to be made up of an incident wave in the dominant mode 
within the coaxial line, with magnetic component B,,(r, z) and a scattered field with 
magnetic components B,,(r, z) for z > 0, and B,_(r, z) for z < 0. Bs, is defined by the 


equation, 


B,, = exp (—itkz)/r for a>r>b, 2 <6, (6) 


With k, > 0, the radiation condition at infinity is satisfied if we take B,_ to be of back- 
ward-travelling wave form within the line, and B,, to have a wave form travelling away 
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from the open end of the line. Thus, in the half-space z > 0 the behavior of By, is given 
by the relation 
RB,. ~ f(g) exp (—tkR) as R- &, 0<¢< 7/2, (7) 


where R = (r° + 2’)'”, g = tan’ (r/z), and f(y) is a bounded function independent of 
Rk. At the edger = a, z = 0, the condition on the total field which ensures the integrability 


of the edge current and hence the uniqueness of the solution [2, 5] is that B, = con- 
stant — O(c°*), and dB,/dz = O(c *”*) as the distance o of the point of observation to 


the edge tends to zero. The problem is now completely and uniquely specified. 
II. To find the field in the half-space z > 0, we define a Fourier cosine transform 
for all r > b by the equation, 
F(p,r) = | B, .(r, 2) cos pz dz. (8) 
The absolute value of this integral for large values of z is found from Eq. 7 to contain 
the factor exp — k; (r° + 2°)'”’, and this factor ensures the absolute convergence of the 


integral. The inverse of Eq. 8 is given by 


B,.(7,2) = 2 | F(p,r) cos pz dp/r, (9) 
for all r > b. 
If we multiply Eq. 2 by cos pz and integrate with respect to z in the range z > 0, 
we find that 
ao ld ; Pi 0B 
(: 7+: + K* — 3)\Kop,7) = (2 
or r or r OZ / <0 
so that for r > a, from Eq. 5, F(p, r) satisfies the equation 
a 1 a “— 
(< s+- : + k* — 1)F, r) = 0, (10) 
\or r or r 
and fora > r > b, writing L(r) for (0B,/dz),-0 , F(p, r) satisfies the equation 
las 
oO lo _ := 
(S s+-~-+K - FO, 7) = L(r). (11) 
\or r or Fr 
Here K = (k° — p’)'’*; that branch of the square root is taken for which when k, = 0, 


K is real and positive when p is real and | p| < k, and for which 7K is real and positive 
when p is real and | p| > k. 

From Ey. 3, d[rF(p, r)]/dr must vanish at r = b. On the surface r = a, B, and dB,/dz 
are functions of integrable square for all z > 0: the continuity of the field components 
B, and E, across this surface therefore ensures the continuity of the functions F(p, r) 


and 6[7F (p, r)]/ér at r = a. From Eqs. 7 and 8 it follows that as r — © 


Fip,r) ~ f(g) exp (—72kR) cos pz dz/R. 
I J I 


For large values of r this integral may be evaluated by the method of steepest descents: 
the calculation shows that 
F(p,r) ~ Wp) exp (—tkr)/r'”” 


where ¥/(p) is a function of p only. 
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The solution of Eq. 10 which shows this behavior for large values of r is a multiple 
of H,(Kr); H,(Kr) is that solution of Bessel’s equation for which as r — © for real or 
imaginary values of K 

(Kr)'?H,(Kr) ~ exp [—i(Kr — 32/4) ][2/r]'”. 
When XK is real this is the Hankel function of H{* (Kr). Similarly, H,(Kr) is a zero order 
Bessel function which for real values of K is to be identified with the Hankel function 


Hy,’ (Kr). 
The solutions of Eqs. 10 and 11 which satisfy the required boundary and continuity 
conditions for F(p, r) and d[rF (p, r)]/dr are, for r > a. 


F(p,r) = +H,(Kr) [ pL(p)[J (Kp) ¥( Kb) — Y,(Kp)J.(Kb)] dp/2H,(Kb), (12) 


b 
and forb <r <a, 


2F(p, r)[Jo(Ka)¥(Kb) — ¥o(Ka)Jo(Kb)]/x = [J\(Kn Y(Kb) — ¥(Ka)J(Kb)] 
[pp H(Ka)(J.(K)¥(K0) — ¥.(Kp)Jo(K0)] dp/Ha(Kb) 
~ [J.(K) ¥A(Kb) — ¥(KnJ(Kb)] on 
-[ pLp[Ju(Kp) ¥(Ka) — ¥\(Kp)Jo(Ka)] dp 
~ [J\(Kr) ¥o(Ka) — Y\(Kr)J(Ka)] 


-[ pLip)[JKp) ¥o(Kb) — ¥(Kp)Jo(Kb)] dp. 


These equations determine the cosine transforms of the magnetic field in the half-space 
in terms of the function L(r) which is proportional to the radial electric field in the 
gapz=0,4a¢>r > DO. ‘ 

From Eqs. 9 and 13 we may write down the total magnetic field in the gap. This is 
Beir, ) = [ aph PERAVAKY — Y(KAIAKO)\ HoKa) 
ic Jo "J (Ka) ¥(Kb) — Yo(Ka)J.(Kb)) Ho(Kb) 

s yee ne . J (Kr) Y¥,(Ka) — Y,(Kr)J,(Ka) 
Lp) [J (Kp) ¥o(Kb) — J J (Kb)] dp -— SSS SS = 
| ple) WS (Ko) Folk (Kp) Jol K)] do — 7 Ka) ¥,(Kb) — Yo(Ka)Jo(K6) 





(14) 
é tate ao . J (Kn) Y (Kb) — Y,(KnJo(Kb) 

° dha ) ‘ ? { b — vo ) oe SS fie sae oe ee — - & sah 

| pLAp) (I (Kp) Fo(KO) — ¥i(Ko)Jo(K0)] do — 7 ey by — ¥(Ka)J (Kd) 


| pL(p)[Ji(Kp) ¥o(Ka) — ¥\(Kp)Jo(Ka)] dp. 


III. In the region z < 0,a > r > b, the field can be described in terms of a set 
of axially symmetric transverse magnetic modes. These modes are associated with the 
set of functions r'”’ ¢,(r), orthogonal in a > r > b, where 


gor) = 1/r, 
and 


g(r) = J,(K,r) Y.(K,a) — Y,(K,r)J.(K,a), 
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for integer values of n > 1. When n > 1, K = K, is the nth positive zero in order of 
magnitude of the function J,(Ka)Y,.(Kb) — Y.(Ka)J.(Kb), and Ky = 0. Foralln > 1, 
(a — b)K, ~ nz [1]. If we put p, to be that value of p which corresponds to the value 
K = K,,, then the magnetic field component of any axially symmetric TM field within 
the line is a linear combination of modes of the form ¢,(r) exp (= tp,z). 

We have assumed that an incident dominant mode B,, of unit amplitude travels 
towards the open end of the line. If R, is the reflexion coefficient associated with the 
nth mode reflected at the plane z = 0, for an incident dominant mode of unit amplitude, 
then we may write the reflected field in the form 


Bs. = >. R.¢,(r) exp (tp,2), (15) 
so that 


6B,_/éz = > ip,Rig,(r) exp (tp,z). (16) 


By matching both the total magnetic field and its normal derivative across the open 


end of the line we find that 


B,.(r, 0) = >> (R, + bo¢,(7), (17) 
and 

Lr) = 2. ip,(R, — bon)¢,(7r), (18) 
where 6,,. = 0 form ¥ n, and 6,,,, = 1 for all integer values of m. Both the functions 


r'?*L(r) and r'”’B,, (r, 0) are functions of integrable square in a > r > b, so that they 
may be expanded in orthogonal series of the functions ¢,(r). 
Thus, ifina >r>b 
Lr) = >> ag,(r), 
n=O 


then 


aa 


| rL(r)y,(r) dr, (19) 


I] 


a,N , 


where 


a 


N, | rly.(r)]? dr. 


II 


|\N, = In (a/b), | 
\. 

N, = 2[1 — [Jo(K,a)/Jo(K,.b)]?]/(rK,)°| 
We can see therefore from Eqs. 18 and 19 that 


a, = tp,(R, — Son); (20) 


and from Eq. 17, that 


. 


| 


6 


rB,,(r, O)e,(r) dr = N,(R, + 4o,). (21) 


Since we are able to write B,, (r, 0) in terms of the coefficients a, , Eq. 21 may be manipu- 
lated to represent an infinite set of equations linear in the a, whose solution will determine 
the field throughout the system. 
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Thus, from Eqs. 14 and 19, assuming that the required changes in the orders of 
integration and summation are permissible, it follows that 








a rB, .(r, O¢,.(r) dr = b a (22) 

b n=0 

where 

’ ait H,(Ka)¢,,(ag,(a) 
C.,=k [ apd Ka! _— = =~; ae 
Jo H,(Kb)\(K° — K,)(K — K,,) (23) 
ees " Q2N,, 5 

. Gi c= r r e h — oe 
[J.(Ka) Y,(Kb) Y,.(Ka)J,(Kb)} 7K? — +K° 


It follows from Eqs. 21 and 22 that the infinite set of equations relating the coefficients 
a, 18 


2 SomNin + AN m/(ipm) = D> CmnQn/k (24) 


for all integer n > 0. 

IV. In the complex p-plane the integrand defined in Eq. 23 has for singularities 
only branch points at p = + k, . It is clear that the path of integration for the integral 
C..» Must be deformed into an equivalent contour: this is necessary in order to avoid 
passing through the branch point p = + k in the limiting case when k; = 0. Thus, in 
the neighborhood of this branch point we can deform the path of integration into a 
semi-circular are of small radius 6 in the upper half p-plane. To examine the contribution 
to the integral C,,,, from this arc we put p = k + 6 exp ig, 0 < @ < 7, where we can 
now take k, to be zero. Since as 6 — 0 the corresponding form for K is 


K = (2k 8)'” exp 7(@ — x)/2[1 + 0(8)], 
we can show that 
a{Jo(Ka)¥.(Kb) — Y,(Ka)J,.(Kb)]/2 = In b/a + Ofé(a — b)]; 


it follows that for m = n = 0, when N, = In a/b, the absolute value of the integrand 
for small values of 6 is 0(1), and the contribution to the integral Coo is 0(6). For m orn 
not zero the contribution may be shown to be even smaller in magnitude. We may 
therefore write C,,,, in the form of a line integral, this being the limiting case for 6 = 0. 
The question of the convergence of this integral is considered in the Appendix. It is 
shown there that C,,,, is uniformly convergent with respect to the parameters a and b 
in the range a > b > O. Since we have defined the branch of k so that for p 
real, 0 < p < k, kis real, and for p > k, k = — ty where uz is real, it follows that 


. aly, (ae, 
( a= [ adK ger H(Ka)yn(a en(4) a a 
k /0 \ H AK b)[K —_ K,, |(K — i 


-[Jo(Ka) ¥ (Kb) — Yo(Ka)Jo(Kb)] + > Np tow 
af i , K (ua) ¢(a)y,(a) 
: ant a? —— es <r 
+ a I I \" ° Ko(ub) (uw + K;,)(w + K;) 
Nnemn 


[To(ua)K (ub) — Io(ub)Ko(ua)] — Ke + a} 
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where K = (k° — p’)'”, wu = (p” — k’)'”’, and K,(z), Jo(z) are modified Bessel functions 
(6). 


V. The first quantity of physical interest to be considered is the energy density at 
large distances from the antenna. From Eqs. 9 and 12 it may be seen that for r > a, 


B,.(r,2) = | cos pzH , (kr) | Lo)[d (Kr) Y,(Kb) — Y,(Kr)J.(Kb)] dp dp/H,(Kb), 


Jb 


so that as r — ©, since the integrand is an even function of p, we have that 


dp exp (ipz — 1Kr) 


B, .(r, z) ~ const. | A : 
a ‘ (Kr)'”"H,(Kb) or 
(26) 


| plAp)[J (Kr) ¥(Kb) — Y¥,(Kr)Jo(Kb)] dp. 


d 


This integral may be evaluated by the method of steepest descents. By putting 


R = (r’' + 2’)'”, o = tan" r/z with O < ¢ < x/2, we find that the saddle point is at 
p = — kcosg, corresponding to the value k = k sin y. The path of integration is deformed, 


without changing the value of the integral, into the steepest descent path, which at the 
saddle-point makes an angle of 7/4 with the line joining the branch points. The value of 
B,, for large values of R can now be given. This is 


exp (—7ikR) “ae . . : 
B,..~ const -_ 2 aE [ pL(p)[J,(kp sin ¢) Y,(kb sin ¢) 
6 


"kRH,”(kb sin ~) Ji i 
, ‘ ; (27) 
— Y,(kpsin ¢)J (kb sin g)] dp, 
and from Eq. 19 after performing the integration, we find that 
const. sin g exp (—7kR) : : 
B, $$. — [J,(ka sin ¢) Y,(kb sin ¢) 
: kRH,”’(kb sin ¢) Jo ' 3 (92) 


a,,ag,(a) 


— Y,(kb sin ¢)J,(ka sin ¢) 5 ea 
* | > Fein? ¢ —K 

The energy density at large distances from the open end of the line can now be found 

from the complex Poynting vector, and it follows that the power flow P in the radial 


direction at large distances is given by 


r sin’ g{ J,(Kasin ¢)Y,(Kb sin ¢) — Y,(Kb sin ¢)J,(Ka sin ¢)}> 
= const. —— a a 3 
R°}[J.(Kb sin ¢)]° + [Y¥o(Kbsin ¢)}°} 29 
Ze } 
j a,,a¢,\a) 
“—~k sin’ o — K,, 


The second quantity which is of interest to us is the admittance of the antenna, 
regarded as a termination of the transmission line. If we go far enough along the line 
for higher order modes to be practically damped and then carry out standing wave 
measurements, we may determine the dominant mode reflexion coefficient R, . The 
terminal admittance Y, associated with this coefficient Ry is given by the equation 


Y, = Y(#, + Dw — Dd, (30) 


where the characteristic admittance of the line Yo = 2rwe,/k In (a/b). From Eqs. 20 
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and 24, it follows that 


Y, = i¥o > Condn/QoNo . (31) 
mom 4 
In an isolated radiating system the only way in which a driving-point admittance may 
be defined is by the ratio of current to voltage at this point. Thus, in this system if we 
define an admittance Y, to be the ratio of the current at the foot of the antenna to the 
voltage applied between the ground and the foot, then it follows, from expressions which 
can easily be found for the current and the voltage, that 


VY, = i¥ob D> D> CamAmen(b)/aoN, « (32) 
m=0 n=0 
The first term of this series is the admittance Y, . 

VI. The exact solution to this problem for general values of the quantities k, a, 
and b involves the solution of an infinite set of linear equations whose coefficients are 
infinite integrals. Even with the help of modern computing machinery the work would be 
exceedingly difficult. We shall first consider, however, the solution for small values of 
the parameter e = k(a — b). If we put o = ka, and p = b/a, then it is clear that « may 
be small if either ¢ — 0, p ¥ 1 or if « ¥ 0, p— 1: the first case corresponds to the physical 
problem of the antenna which is thin in comparison with the excitation wavelength, and 
the second to the problem of the thick antenna with a low impedance line feed. We will 
exclude, however, the case for which p — 1 and « — 0, this being the case of the thin 
antenna with a low impedance feed. 

For the limiting cases ¢ > 0, p ¥ 1, o ¥ 0, p > 1 the magnitudes of the functions 
which appear in the infinite set of equations are derived in the Appendix. To solve the 
problem when ¢ — 0, we first assume that in the equations 24, a, = 0 for n > 1. Then 
from the first equation of the set, it follows that 


a = —2ik, as pl, o¢ ~ 0, 
or o— 0, p¥l. 
If we take this first approximation, we find from the (n + 1)th equation that 


oC. 
- Ri: 
= 0(1) as pol, o £0, or oO, pl. 


Using these results for the orders of magnitude of the a, in the complete set of equations, 
we find that as e —- 0 
a, = —2ik{1 + 0(6], 
a, = 0(1), ry > i. 
We now find in the expressions for the physically interesting quantities that 
ao = — 2k is the only coefficient to make a significant contribution. This value for a> 


is that which would arise in analysis of an open-circuited line ignoring end effects. 
From Eq. 19, the unknown function may now be given in the form 


L(r) = —2ik/r. 








V. M. PAPADOPOULOS [Vol. XVII, No. 4 
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ADMITTANCE Y= G+iBle=k(o-b), p=b/a , Y¥, = 24we,/k In (a/bd)) . 

| 

<j G/Y (e€ =0.1) 
| 

4 

3 = 
| 

+ 
B/Y, (€ = 0.1) 

| t 
ho. B/Y,(¢ 0.05) 
\r=-- ef *U 

o| nafs L, 
fe) 0.2 0.4 0.6 0.8 , 


Fic. 2 


The corresponding admittance of the radiating system, from either definition in Eq. 
30 or 31 takes the form 
Y = ikY Coll + O(€)]/No . 


The value of Y in the limit as c — 0, p ¥ 1 is O(e). The limit as p > 1, o ¥ O does not 
exist because the second derivative of Coo with respect to p does not exist. In this case 
the strongest statement to be made is that the normalized admittance Y/Y, is zero. 

Numerical values of the admittance for small values of ¢ are given in the accompany- 
ing graph (Fig. 2). The normalized admittance approaches unity as p — 0 for fixed small 
values of c. The open-circuit approximation used to calculate results is not therefore ap- 
plicable for small values of p. From Equation 29, when « — 0 we find that 





2a [J (osin ¢)Y(epsin ¢) — Y,(o sin ¢)J,(cpsin ¢)]’ 
P mw oe, SLE eee eee [1 + 0(6)]. 
uk’r [Jo(op sin g)|} + [Yo(epsin ¢)] 


For r > a at large distances from the end of the line, where tan gy = r/z this expression 


may be simplified in certain cases. 


In the case when co > 0 p ¥ 


P ~& [In pl’, ‘uk’r? { [In (op sin g)/2]’ + 2°/4}. 
In the case when p —> 1, o * 0 is more interesting. Thus 
P %& 801 — p)’/pk’r’x’{[Jo(op sin ¢)]’ + [Yo(opsin ¢)]}’}, 
and when sin g  r/z for small enough values, r > a, 


P ~& 8w(1 — p)?/pk’r’e’{4[In (or/2z)/r]? + 1}. 


This expression shows that we may expect a behavior for P like [1/r In r]’ for short wave- 
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lengths when ¢ = ka = 0(z) for large values of z. More generally, if, when p — 1 and sin 
v is not small, we take values of o which are large, but not so large that the relation 
a(1 — p) = e— O is not satisfied, we find that 


P & 4dwe(1 — p) sin ¢/urk’r’. 
VII. In the general case, for « > 0, we shall prove two results. The first is 


Theorem 1. The infinite set of equations given in Eq. 24 has for a solution at most one 
set of values a, tf « = ka(ab) > 0. 


It will be recalled that a, is a coefficient in the orthogonal expansion ina > r > b 
of the function 7'” L(r). 


Suppose that we have fixed rr < « < x (r + 1) so that p,, (0 < m < r) is real and 


iDm = Om (m > r+ 1) is also real. The hcmogeneous set of equations corresponding to 
Eq. 24 is 


KN ,.Q,/1)m = > Cim1m , for m> 0. (33) 
Putting a,, + 78, = an, Cnn = Kan + tlm, We find that 
ie (KmnOQin — LinnBan) = kKNpBn/Pm , for m <r, 
= kN,an/Qn, for m>r+1. 
Dx (KumBn + LimnOQn) = —KN,Qn/Pm, for m<r 


= kN .8./Ga ; for m>r+l. 


It follows that 


—Bn >. LnnBn — Om >, Lat, = 0, for m>r+1, (34) 
= kN,(a2, + B.)/pn , for m<r, 
so that 
k > Nala, + Ba)/Pm = — DD (BalimnBn + mL mr Qn) + 
Now 
] 9}. [ ¢m(a)en(a)K*a*[Jo(Ka) ¥ (Kb) — Yo(Ka)J(Kb)) 5 
ao (K — K2)(K — K){[J(Kb] + [YoKbF} “”’ 


so that the result of changing the orders of integration and summation in Eq. 34 is 


: . | ’ dp K*a° 
f - an f ) = —2 we _ 
2 N alae + Bn)/Pm 7 [ {[J(Ka)” + Y(Kb)}} 
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The right-hand side of the equation is < 0. The left-hand side of this equation is, how- 
ever, a non-negative quantity so that each side of the equation must be zero. It follows 


that 
Qn ,Bm = 0 forall m>O. 
Given the infinite set of equations, it is necessary to know whether an approximate 


solution can be found by solving only the first N + 1 equations for the unknowns 
(n < N), taking all the remaining unknowns to be zero. We shall therefore prove the 


following 


Theorem 2. (i) The determinant of the finite set of equations 
N 
kN, (2 Som + an /tDm) = tk Canis ’ 0 < m < N, (35) 


does not vanish for any positive value of «€ 
(ii) The unique solution of the finite number of equations is the set {aX},O <n < N, 


and lim (a, — a“) = OasN > o. 
The first part of the theorem is proved in the same way as in Theorem 1. For the 
second part we find from Eqs. 24 and 35 that 
r. N , 
KNn(Qm — On)/Pm = Dy Cum(Qn — A) + DO Camden 
n=0 + 
Now from 19, using Parseval’ Theorem, since r'”* L(r) is a function of integrable square 
ina > r > b, we see that 
| riL@P dr = | a, |? N, = 01), for «> 0. 
vb 0 
Since VN, = O(r-’) as n — ©, a necessary condition for the convergence of this series 
is that asn — o, for 6 > 0, 
| a, | = O(n'~’). 
It is shown in the Appendix that the infinite integral C,,,, is absolutely convergent for all 
e>0.Asn— @- |C,,, | = O(n™*) so that 


oe =] YS Canta | = O(N) 
n=N l i 


if N is large enough, for all m > 0. 
Thus 


kN (Gm — ay) . 
alts — 5) _ 5° C_(a, — a!) | = of, (36) 
1D m pe 
and since the determinant of the finite set of equations is non-zero, the solution of 36 
must be of the form 
N oe 
a, — a, = Oa,,). 


Therefore, we may approach as close as we like to a solution of the infinite set of equations 
if we take N large enough, and from Theorem 1 this solution is a unique one. 
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Theorem 2 shows that we may obtain a useful approximation to the solution by 
solving for the first N + 1 unknowns, and we are able to use a variational method to 
hasten the convergence to a correct result. This approach may be made in a more general 
problem in which there may be set up any number of propagating modes within the 


coaxial line. 


APPENDIX 
We transform the integrand in Eq. 23 by putting p = k cos 6. To correspond to the 
branch of K, we must take the path of integration in the complex @-plane as shown in 
the diagram (Fig. 3). 
We are principally interested in the behavior of C,,, for small values of the quantiy 
e = k(a — b). We have to consider separately the possibility of the ratio p = b/a approach- 


ing unity with = ka not small, and also the case in which o is small and p is not close 


to unity. 
From the definition of the functions ¢,(a) in Sec. III, it follows that as e — 0 


ag,(a) = —2e/nkr*, for n> 1 
= }, for n = 0. 


The functions N,, defined in the same section behave in the following manner for 
small values of e. Thus, for any value of ¢ 
( i ~~ I, 
No. =hp= 66 if 

Oe) if pl, 

and 
N, = 2{1 — [Jo(K,a)/Jo(K.b)P}/(rK,)’, 
Qe” 


eo. 
7 wk n? ' 


— [Jo(nx/1 — p)/Jo(nmp/1 — p)}}, 


=O) if p¥1 «0, 
=O) if pol o #0. 
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We shall first consider the convergence of the integral C,,,,. On the path C we have 


already established the bounded nature of the integrand near the origin, since we have 
examined its behavior in the p-plane near the branch point p = + k. This property 
is independent of the values of the parameters o and p. On the negative imaginary axis 
the integrand J,,, is a real function involving modified Bessel functions of real argument. 

We must examine the uniformity of convergence of the integral {% J,,,d¢ with respect 
to the parameters o and p. Let us take a small positive number a» as close as we like to 
zero. Then, for 0 < p < 1, ¢ > a > 0, we may choose g > sinh™ (1/o,), and for 
¢, > ¢ we may use the asymptotic expansion of the modified Bessel functions. Since 


for all n > 0, K,, > 0, it follows that 


is) 


ag,\ajag,,\ a) 


| te |< | om E Oun Tt 7 - —(1 + O(esinh g) » | 
1 Je J., sinh ¢ 2o sinh ¢ 


a 


4 


Since for 6 > ¢; , ¢ om, 8 > 1, it follows from the convergence of the integral {?, cosech 6d@ 
that C,,, is unformly convergent with respect to ¢ and a in the given range. This is true 
for all m, n > 0. 

Now consider C,, . The integrand vanishes for c = 0, p ¥ 1 ora ¥ 0, p = 1. It is 
clear that for p ¥ 1, a isa factor of Coo , and for a = 0, (1 — p) is a factor. The derivative 
of the integrand with respect to p is uniformly integrable in the same range of parameters, 
and this vanishes also for p = 1; it follows that for ¢ ¥ 0, pl, Coo = O(e€’) while for 
p © 1,¢—0Co = O(e). The second derivative with respect to p does not exist if p = 1, 
o ~ 0, nor does the first derivative with respect to o when o = 0, p ¥ 1. 

The magnitude of C,,, for at least one of m or n not zero is determined by the magni- 
tude of the integrand on the path C at a finite distance from the origin. Since K,, ~ nrk/e, 
we find that 

Con = Cu = O(c’) for p— 1], o ~ 0, 


= O(c) for p¥1, ¢—0, 
and that for both m and n not zero, 


C 


ll 


O(e’) for p—1, o ~ 0, 


mn 


ll 


O(c‘) for o—0, pl. 
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—NOTES— 


ON THE WORK OF A FORCE CROSSING A BEAM * 
By L. MAUNDER (University of Edinburgh) 


1. Introduction. Timoshenko [1] has discussed the paradox that during its transit 
across an elastic beam whose ends are supported at the same horizontal level, as shown 
in Fig. 1, a vertical force apparently performs no net work, yet leaves the beam in a 
state of oscillation. His explanation, however, refers to a different but closely related 
problem, in which the force is always directed at right-angles to the deflected beam. 
A horizontal force component is thereby introduced which evidently is capable of doing 


work. 














Fic. 1 


Lee [2] shows that this work equals the sum of the strain and kinetic energies possessed 
by the beam on departure of the force. But, pointing out that this result does not explain 
the original difficulty, he returns to the case of the purely vertical force and shows that, 
contrary to appearances, the work which it does on the beam is not zero, but equals the 
energy retained by the beam. This result then resolves, in principle, the original paradox. 

In view of the usual interpretation of the work of a force, however, it remains curious 
that the purely vertical force does positive work during transit between points on the 
same horizontal level. It therefore seemed worthwhile to consider the problem further, 
not only because of the physical mechanisms involved, but also because of its bearing 
on the definition of work. 

2. Work relations. As in [2], we define the rate at which a force does work on a 
body as the scalar product of the force vector and the velocity vector of the particle of 
the body at the point of application of the force. This definition permits association of 
the work done on the body with the change in its kinetic, strain and potential energy. 
But it does not always describe fully the transmission of work from one physical system 
to another, as may be illustrated by the system of Fig. 2. 








Fa Ad Ma F 
B ----2V, 
77774777 7 
Fig. 2 





*Received July 30, 1958; revised manuscript received October 27, 1958. 








L. MAUNDER [Vol. XVII, No. 4 


Block B slides on a frictionless horizontal plane with velocity Vs, : the center of the 
circular disc moves horizontally only and the disc rotates such that the velocity of 
particle A in contact with the block is V4, . As the normal components of the contact 
forces between disc and block do no work, only the frictional components fF, and F', 
are shown. From the preceding definition we find that /, does work on the block at the 
rate (— F,V3,) and that F,4 does work on the disc at the rate (F,V 4). This description 
can usefully be expanded, however, by stating that FP, extracts work from the block at 
the rate (/';V,), of which part is transmitted to the disc at the rate (F,V,4) and part is 
transformed into heat by the force pair F(F, = Fs, = F) at the rate F(V, — V,). 

If V, is zero, when the figure might represent metal-cutting with a fixed tool at A, 
work is extracted by the contact force from the traveling work-piece B at the rate 
(FV). All of the extracted work is transformed by the contact force pair and no work 
so that the disc rolls on the block, work is 


is done on the fixed tool. If V, = Vz = V, 
extracted by the contact force from one body at the rate (FV) and is transferred entirely 
to the other. In this case, no work is transformed into heat. 

In general the flow of work from one system to another should be studied through 
the action of the force pair or pairs exerted between the systems. When tangential 
forces occur at moving surfaces, the foregoing definition of work can lead to the result 
that the work extracted from one system differs from that transmitted to the other. In 
certain cases, the difference is transformed into energy associated with microscopic 
phenomena, such as thermal energy. 

3. The physical system. The application of a purely vertical force to the deflected 
beam requires that the surface of the beam be capable of sustaining tangential forces. 
They may be associated with sliding or rolling. We consider both possibilities. 

A particle of negligible mass can slide across the beam in the prescribed way under 
the action of an applied constant vertical force and the equal and opposite reaction of 
the beam, assuming that the frictional properties of the surfaces allow the latter to be 
developed. The total work done on the particle by the applied force is expressed by the 
integral 

‘ f ' (2 m y 2H) dt (1) 

0 ot C2] cave 


in which the integrand is the vertical component of the velocity of the particle. As the 
net vertical displacement of the particle is zero, the work (1) is zero. Since the particle 
stores no work, the total work extracted from it by the contact force applied by the 
beam must also be zero. Now the (zero) extracted work is partly transferred to the 
beam and partly transformed into heat by the frictional components of the contact 
force pair. The net work done on the beam is expressed [2] by the integral 


al/V e 
VY 
P | (24) it 2 
0 ot z=Ve , \ ) 


which, being equal to the energy retained by the beam, is necessarily positive. Hence 
the work transformed into heat, namely, 
“¥ fay 
Pv [ (24) dt (3) 
0 Ox z=Vt 


< 


must be negative, from which we conclude that the assumed frictional characteristics 
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cannot exist even in principle. The detailed explanation is that the assumptions entail 
negative friction, i.e. friction tending to increase the relative velocity between particle 
and beam, when (dy/dr),-y, is negative. Sliding friction must therefore be discarded 
as a possibility. 

The alternative is rolling contact, of which the simplest case is the rolling of a circular 
dise of negligible mass and radius r. Its point of contact with the beam will have the 
prescribed horizontal velocity V if the velocity of its center has horizontal and vertical 
components V and [(dy/dt) + V(dy/dr)],.y. respectively, and if the dise has a clock- 
wise angular velocity 


(V/r)[1 + (dy/dx)7 7), . 


This motion requires that certain forces be applied to the disc in addition to the pre- 
scribed vertical reaction P of the beam. They are equivalent to a vertical force P applied 
at the center of the dise together with a counterclockwise couple of magnitude 
Pr{(dy/dx)/[1 + (dy/dx)*]'7} eve . 

As the dise stores no work, the work done on it by the applied central force and the 
applied couple must be extracted fully by the contact force of the beam. Moreover, 
the work so extracted is transferred entirely to the beam, since no work is transformed 
at a rolling contact, as noted previously. Now the applied central force and the applied 
couple do work on the dise at the rates 


P{(dy/dt) + V(dy/dx)|,-y. and —PV(dy/dz),-y: (4) 


respectively, as follows from definition and the prescribed motion of the disc. Hence 
the rate at which work is transferred by the contact force to the beam is the sum of the 
rates (4), namely, 


P(dy/dt),-ve » (5) 


If we regard the central vertical force and the couple as external forces which are 
applied to the system of dise and beam, it is noteworthy that the total work done on the 
system by the central force is zero. The work done on the system by the couple equals 
the work done on the beam by the contact force, i.e. 


oy fo 2u) —- (24) 
—Py [ (22 ae di =I i at) .-v: dt. (6) 


4. Conclusion. A full description of the work of a force requires specification of 
the physical nature of the systems concerned. A physical inconsistency appears if the 
vertical force is considered to be applied to the beam through a sliding body, but it 
can be applied through a rolling body. In the simplest case of a rolling circular disc of 
negligible mass, the external forces which must be applied to the disc are equivalent 
to a constant vertical force applied at the center of the disc and a couple. The vertical 
force does no net work on the disc during transit of the beam, and the couple does an 
amount of work on the dise equal to that retained by the beam. 
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AN IMPLICIT, NUMERICAL METHOD FOR SOLVING THE 
n-DIMENSIONAL HEAT EQUATION* 


By GEORGE A. BAKER, JR. (Los Alamos Scientific Laboratory, Los Alamos, New Mexico) 


1. Introduction. The work of this paper represents an extension of the work of 
Baker and Oliphant [1] to the n-dimensional, linear, heat-flow problem with arbitrary 
spacing between the mesh points. We give an implicit scheme which we are able to solve 
exactly, and we prove that it is unconditionally stable. We treat only the linear case 
because the generalization to the non-linear case may be treated in precisely the same 
way as was done by Baker and Oliphant and the same conclusions hold in the 
n-dimensional case concerning this generalization as in the two-dimensional case. 

2. The difference approximation. The basic partial differential equation we wish to 
consider is 

5 , 2 , 00 
Vo= Sx,) +A 0+ B af? (2.1) 


where V’ is the n-dimensional Laplacian with respect to x, x is an n-dimensional position 
vector, ¢ is time, and S is a source function. Not both A and B may be zero. By setting 
B = 0(A # 0) we may solve at one step for the asymptotic solution of (2.1). We wish 
to approximate this differential equation by an implicit, difference scheme to allow 
approximate numerical calculation of the function 6. We shall define our mesh by n 











sequences 
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which are the values of the coordinates of the mesh points in an n-dimensional vector 
space. 
Let us define the matrices 
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where 6) is the Kronecker delta. Let us approximate 
06 oe : 
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C 
and define 
B = —[A’ + 3B’*/(2 Ad)]. (2.5) 


We may now consider the nth order tensor (suppressing the various Kronecker 6’s) 
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It is evident that 
2 


>> D'',6;, represents a (2.7) 
Ik OX, 


to within terms of the second order in Az, . Hence (2.6), neglecting 4th and higher order 
derivatives, represents to within second order in Az, the left hand side of 


86+ 770 = S(x, ) + A?0+ B’ - + 86. (2.8) 


By definition 8 has been chosen so that the right-hand side of (2.8) is independent of 
6(t) and depends only on @(¢ — At), 6(¢ — 2At), and a known function, S. We can solve 
the implicit equations (2.6) for 6;,...1,...7, , a8 86,* + Dj* is a triple diagonal matrix. 
We can factor it as 


Bs, + DF, = > witb, , (2.9) 
Le 
where w is a lower triangular matrix. The matrix elements are given by the recursion 
relations 
wii =0 if i, —-L, 0,1 
bo =0 if I, —L, ¥0,1 


bt = 1 
(2.10) 
Je-1 \Je-l 

wy, = D3; 

wy, = B+ Dz; — ws'b3t-, (bo = 0) 

b3i*? = Dzi*'/wst . 


Let us now define 


Gr, cceLygeoeL, = = (17 a ° (2.11) 
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Te 


Then, using the difference approximations (2.4) and (2.6), we may represent (2.1) by 


po (TI | ee ee (2.12) 
where the a’s are known quantities. Due to the triangular nature of the w’s we may 
proceed from low index numbers to high ones and solve (2.12) for the g’s by straight- 
forward elimination. Then due to the triangular nature of the b’s we may proceed from 
high index number to low ones, and solve (2.25) for the 6’s by straightforward elimination. 

3. Stability. We shall proceed with a proof of the unconditional stability of the 
method described in Sec. 2 by a method closely related to those employed by Baker and 
Oliphant [1]. We consider only the problem with S independent of ¢. As Eq. (2.1) is 
linear it will be sufficient to consider the solutions of the difference equation analogue of 


770 = A’6 + B - (3.1) 


with homogeneous boundary conditions. This equation results from the subtraction of 
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the asymptotic solution from the difference equation analogue of (2.1). Let us seek a 
solution of the difference equation analogue of (3.1), such that, for some Z, , 


Ov, x, t) = wv, x)Z/™'. (3. 


- 
bo 
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Equation (3.1) as expressed by (2.4), (2.5), and (2.6) is for a component (3.2) 
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and an operator A, , such that 
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One can now show that the solution of (3.3) is the same problem as finding the extrema 
of F + dG subject to (3.7), and homogeneous boundary conditions. As we have 


IT M, = M (3.8) 


mesh-points, the theory of quadratic forms [2] assures us that since 8 is negative, there 
exists a complete set of m orthogonal vectors (67,...7,...;,(¥)) Which satisfy (3.6) and 
(3.7) and therefore (3.3). As \ can be shown to be the negative of F, we must have 


(3 — 4Z,' + Z,*)/3 = At(\ — 2A°/3)/B’ < 0. (3.9) 
It easily follows [1] from (3.9) that for both roots 
ae ae ie (3.10) 


Hence, the amplitude of all the 6(v) solutions decays to zero as time goes to infinity. 
Since the 6(v) form a complete orthonormal set, any initial value problem can be expanded 
in terms of them, and hence we may conclude on the basis of (3.10) that the method is 


unconditionally stable. 
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HELICAL FLUID FLOWS* 
By ROBERT H. WASSERMAN (Michigan State University) 


Introduction. Potential helical flows have been completely described by G. Hamel 
[1]. Nemenyi and Prim, and N. Coburn have obtained some Beltrami helical flows 
[2, 3]. A simple description will be given here of all steady incompressible helical flows 
and all steady compressible helical flows with entropy constant along stream lines. 

The equations of helical flow. The class of compressible flows to be considered 
here are those governed by the following differential equations in which o* is the velocity, 
p is the pressure, p is the density, S is the entropy, and é* is the unit tangent vector 


along the stream lines; 


VY -pv* = 0 (continuity equation) (1) 
(pv*- 7 )vo* = — Vp (equation of motion) (2) 
p = p(p, S) (equation of state) (3) 

i*-7S = 0 (entropy is constant along stream lines) (4) 


For incompressible flows we must satisfy Eqs. (1) and (2) and the special case p = 
constant of Eq. (3). 

Now we assume that the flows are helical; i.e., the stream lines of the flow are parallel 
helices on coaxial circular cylinders. Such flows have the property VY -¢* = 0. This may 
be seen immediately if we introduce cylindrical coordinates r, 6, 2 and decompose ¢* 
according to 


t* = sin B6* + cos Be*. 


Note that the angle @ of the helices is in general a function of r. 
With the condition VY -t* = 0 the continuity equation reduces to 


i*-Y7 In pq = 0, (5) 


where q is the magnitude of the velocity. Also, from Eqs. (3) and (4) 


) 
i*-Tp — Pa i*-Tp = 0. (6) 


Finally, we write Eq. (2) in intrinsic form [4]: 
Pipa eek, G 
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n*-Tp = —pq'k, (8) 


b*-Vp = 0, (9) 


where n* is the unit principal normal vector of the stream lines, b* is the unit binormal 
vector of the stream lines, and x is the curvature of the stream lines. Then Eqs. (6), 
(7), and (8) give us the result that either 


(a i*-Vg = *-Vp = t*-Vp = 0 
or 
‘ 2 9 A(S) 
(D) q = Me and p= B(S) — As . 
dp p 


where the functions A(S) and B(S) are restricted only by the condition A(S) > 0. 

Clearly, the alternative (b), in which the equation of state includes that used by 

Chaplygin, and Karmd4n and Tsien, only occurs in the compressible case. 

The fact that the helices are geodesics on the cylinders r = const. means that n* 

is normal to these cylinders. With this, and using Euler’s equation for the normal curva- 
ture of a curve to evaluate x, Eq. (8) becomes 

op at 2 Sin’ Pp (10) 


( 
ar PY , 


The two classes of helical flows. In case (a) Eq. (10) reduces to 


dp > sin’ B 
— <— - 
dy Y 


If in this equation it is understood that p and q are constant along stream lines (i.e., 
condition (a) holds), then it embodies all the conditions for helical flows in case (a). 
That is, there is a helical flow corresponding to any set of functions p, p, g, 8 satisfying 


this equation. 
In case (b) if we introduce the coordinates 


a=2z-+ @rtan 8B, 
y =z — Or cot B, 
r=Tr, 
then Eq. (10) may be written 
Op i d Op s eee 
P< + (a — y) sin B cos B (r tan B) ~ = (B —- p) sin’ Bp. 
or dy 0a 


If in this equation it is understood that p is only a function of r and a, and B is only 
a function of r and y then any solution gives a helical flow. As an example of a family 
of solutions we have those in which 

B const. 


rtan B = A = const. 


(rn? + 9°*)** _ 
ann amend i”. S 


7 


p=B- 


where F is an arbitrary positive function of a. 
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Some general properties of flows of case (a). In case (a) which applies to all fluids 
except those having the special form of the equation of state of case (b) our basic equa- 
tion is very simple and we can readily obtain some general properties of helical flows 
from it. 

Thus we see that 

I. p is a non-decreasing function of r, and it is easy to construct examples for which 
it is either bounded or unbounded 

II. p, g, and 8 may, in general, either increase or decrease with r. Moreover, p and q 
may vary from stream line to stream line on a cylinder as long as pq’ is constant on 
cylinders. 

III. In certain important special cases the variation of p, g, and @ is restricted. Thus, 
for example, in the incompressible case if the Bernoulli function is constant then g must 
decrease with r. In the compressible case, if the entropy and the stagnation enthalpy 
are constant, then g must decrease with r 

IV. For a polytropie gas our basic equation can be written quite simply in terms of 
the Mach number. We find that the Mach number may in general either increase or 
decrease with r. However, if the entropy and the stagnation enthalpy are constant then 
the Mach number must decrease with r. 

V. With respect to the vorticity we note in particular 

(i) in contrast to the other quantities of the flow [in case (a)] the vorticity can 
vary along the stream lines. 

(ii) when the vorticity vector is normal to the stream lines the only possible stream 
line pattern is given by cot 8/r = const. which is that of the simplest helical 
flow; that obtained by normal superposition of a potential vortex and a uniform 
rectilinear flow. 
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A NOTE ON THE PAPER OF MILLER, BERNSTEIN AND BLUMENSON* 
By F. A. J. FORD (Admiralty Research Laboratory, Teddington, England) 


In connection with the above authors’ recent paper [1] on generalized Rayleigh 
distributions, it may be of interest to show how the cumulative distribution function 
of the first order distribution of such a process can be expressed in terms of tabulated 
functions. The cumulative distribution function in question is 


»R »\ NV/2 
["4 (2) exp [-(R* + 4% Wolew-nva(“4) dR, 


which may be written as 
(2y.)"7_ f R A \ 
A Pn\(2y,)' -* Gay) 
where 
p(x, y) = | 2u"*' exp [—(uw? + y’)]I,(2uy) du n = (N — 2)/2. 


When JN is odd the modified Bessel function J, can, as is pointed out in [1], be ex- 
pressed in terms of Hyperbolic functions, so that p,(x, y) can be expressed in terms of 
Error functions. On the other hand, if N is even (so that n is integral), p,(z, y) can still 
be expressed in terms of tabulated functions. 

For, on integration by parts, it follows that 


p(x, Y) = Ypn-i(x, y) — x” exp [—(a° + y) J, (2zy) (n > 1). 
Thus 
2. [y"p, (x,y) —y ‘Dn vin a y)| = D(r,y) — y"polx, y) 
n—-1 
= —exp [—(x + y)] p x” "y'I,..(2ry). 
This result may be written 
p(x, y) = y"p(x, y) — exp [—(a? + y’)] > 2’y" "I, (2zy). (1) 
1 


Now p,(x, y) has already been tabulated (a list of tables is given in [2], so that by 
using Eq. (1) numerical values of p,(z, y) can be conveniently obtained, at least for small 
values of n, from the tables of p,(x, y) and of the modified Bessel functions. 

The author is indebted to the reviewer for suggesting a simpler derivation of Eq. 


(1) than that originally presented. 
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Correction to my paper 


THE RELATION BETWEEN THE FLOW OF NON-NEWTONIAN FLUIDS 
AND TURBULENT NEWTONIAN FLUIDS 


Quarterly of Applied Mathematics, XV, 212-215 (1957) 
By R.S. RIVLIN (Brown University) 


The last sentence should read: 
“Certainly it appears quite likely that a phenomenological theory of the type considered, 
in which the stress in an element of the turbulent fluid (large compared with the eddy 
dimensions) is supposed to depend only on the kinematic variables in that element 
or on the velocity-gradient history of that element, will not be entirely adequate as a 
complete description of the flow properties of the turbulent fluid, since eddies can 
diffuse from one point of the fluid to another.” 
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Analysis of industrial operations. Edited by Edward H. Bowman and Robert B. Fetter. 
Richard D. Irwin, Inc., Homewood, Illinois, 1959. xi + 485 pp. $7.95. 


A collection of 27 papers (published, in the main, within the last four years) describing actual appli- 
cations of quantitative methods to the analysis of industrial operating problems. This material was 
collected for use in a graduate course in the MIT School of Industrial Management. There is a bias 
in the selection of papers towards discussions of the fitting of mathematical models to actual problems, 
rather than the development of new models or new methods. 

The reader is presumed to have a working knowledge of differential calculus and analytical statis- 
tics, and some understanding of linear programming, waiting line theory and the techniques of Monte 
Carlo simulation. 

The papers are divided into five classifications: Applications of linear programming; Other pro- 
gramming applications; Waiting line applications; Applications of incremental analysis; Total cost 


and value models. 
Bruce CHARTRES 


Elementary decision theory. By Herman Chernoff and Lincoln E. Moses. John Wiley & 
Sons, Inc., New York, and Chapman & Hall, Ltd., London, 1959. xv + 364 pp. $7.50. 


The day of determinism is past, even in the field of engineering, and it is rather amusing to remark 
that the determinist has inevitably brought it on himself. By constructing more and more realistic 
models of the physical universe, of greater and greater complexity, he has slowly forced himself into 
@ position where the only escape from stalemate lies in the introduction of new rules into the game. 
The result is that we must coat ignorance of true cause and effect and inability to handle the resultant 
equations, even if we knew them, by the soothing, albeit sophisticated, salve of probability theory. 

The modern engineer who is interested in feedback control processes, in automation, in processes 
involving adaptive control, and in any number of significant uses of digital and analogue computers, 
must acquaint himself with the mathematical theory of decision processes. He must understand the 
difficulties and advantages of uncertainty, and be familiar with the various tricks that are used to 
construct precise mathematical models of imprecise situations. 

Unfortunately, many of the standard works in the field of decision-making start from a high mathe- 
matical level and continue upward. The result is that they are not suited to the needs of the beginner. 

The book by Chernoff and Moses is a book designed specifically for the scientist interested in de- 
cision processes, but possessing only a rudimentary mathematical background. The authors have 
succeeded admirably in what they have set out to do. Carefully and slowly, they illustrate how mathe- 
matical models of decision-making under uncertainty are formulated. They discuss the choice of a 
criterion function, the choice of possible policies, the description by means of state variables, and so on. 

All of this is done in a very pleasant, readable style, with numerous examples, and much important 
discussion. Most important, they are honest with the reader in constantly pointing out pitfalls and 
limitations. Their presentation illustrates very clearly that fundamental ideas can be transmitted with- 
out the pseudo-abstraction that befogs the Bourbaki and their devotees. 

In the second half of the book, they give an introduction to classical statistics. Throughout, their — 
aim is to follow the guiding idea of Wald who conceived of statistics as decision-making under un- 
certainty. In their presentation, they have followed the teaching of the late M. A. Girshick. 

The book is recommended to mathematicians, physicists, engineers, and so forth, who wish to 
obtain a relatively painless introduction to this modern field, either for their own studies, or merely 


to keep abreast of current intellectual activity. 
RIcHARD BELLMAN 


























